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Positron  Emission  Tomography  (PET)  is  a nuclear  tomographic  imaging  tech- 
nique that  measures  metabolic  activity  within  the  human  body.  A new  integral  equa- 
tion that  models  PET  is  presented,  and  it  is  reduced  to  a set  of  infinite-dimensional 
matrix  equations  involving  an  orthogonal  series  expansion  of  the  kernel.  It  is  shown 
that  this  new  model  is  equivalent  to  the  standard  Shepp-Vardi  model,  but  the  kernel 
of  the  integral  equation  takes  a particularly  convenient  form.  In  particular,  it  is  re- 
vealed that  the  kernel  is  superharmonic  and  its  Riesz  Decomposition  is  derived.  This 
provides  insight  into  the  kernel’s  orthogonal  series  expansion  that  reveals  the  matrix 
structure  and  can  be  used  to  compute  the  matrix  elements.  A fast  reconstruction 
algorithm  is  devised  based  on  the  matrix  equations. 
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CHAPTER  1 
INTRODUCTION 

1.1  Introduction  to  PET 

Positron  Emission  Tomography  (PET)  is  a nuclear  tomographic  imaging  tech- 
nique that  measures  metabolic  activity  within  the  human  body.  This  is  in  contrast  to 
magnetic  resonance  imaging  (MRI)  and  X-ray  tomography,  which  provide  primarily 
anatomical  and  morphological  information.  As  Ter-Pogossian  notes,  PET’s  ability 
to  measure  biochemical  activity  has  important  diagnostic  and  research  implications 
since,  “...any  biological  activity  stems  from,  and  is  accompanied  by,  regional  biochem- 
ical changes  in  the  organ  or  structure  in  which  the  biological  activity  takes  place...” 
( [35],  p.l71).  Menard’s  [27]  study  of  regional  brain  function  is  an  example  of  a 
research  application  of  PET.  With  respect  to  pathology,  PET’s  clinical  uses  include 
location  of  tumors,  measurement  of  myocardial  perfusion  and  assessment  of  cancer 
treatment. 

The  fundamental  observation  behind  PET  is  that  compounds  in  a living  sys- 
tem can  be  tracked  by  labeling  them  with  radio-isotopes.  Specifically,  the  distribution 
of  positron  emissions  is  proportional  to  the  nonhomogenous  absorption  of  compounds 
within  an  organ.  This  absorption,  in  turn,  is  a direct  measure  of  biochemical  activ- 
ity. Compounds  such  as  water,  ammonia,  glucose,  hemoglobin  and  palmitate  can 
be  labeled  with  short-lived  isotopes  such  as  oxygen  15,  nitrogen  13,  hydrogen  3 and 
carbon  11.  The  selection  of  a compound  is  determined  by  the  particular  metabolic 
process  of  interest.  For  example,  labeled  glucose  is  normally  used  to  study  brain 
activity. 
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Following  Anderson  [1],  the  steps  in  PET  are: 

1.  Select  a compound  and  label  it  with  a positron-emitting  radio- nuclide. 

2.  Administer  the  labeled  compound  to  the  subject. 

3.  Estimate  the  distribution  of  the  compound  based  on  emission  data. 

4.  Use  the  estimated  distribution  to  determine  parameters  of  physiological  models. 

In  order  to  address  step  3,  which  is  the  focus  of  this  work,  we  must  describe  the 
emission  process  and  detection  mechanism  in  more  detail. 

A positron  that  is  emitted  travels  a short  distance  until  it  collides  with  a nearby 
electron.  The  products  of  the  ensuing  annihilation  are  two  511  KeV  high  energy 
photons,  which  propagate  in  random,  nearly  opposite  directions.  The  angle  between 
them  is  not  quite  180°  because  momentum  must  be  conserved  and  the  net  momentum 
at  the  time  of  collision  is  normally  nonzero.  These  photons  can  be  measured  by  a 
ring  of  detectors  that  encircles  the  organ  of  interest.  An  annihilation  is  known  to 
have  occurred  in  the  tubular  area  defined  by  two  detectors  when  they  register  nearly 
simultaneous  detections.  This  scheme  is  illustrated  in  Figure  1.1.  The  point  marked 
“x”  represents  the  emission  point.  The  line  through  x represents  the  propagation  of 
the  photons.  When  the  two  detectors  register  nearly  simultaneous  counts,  it  signals 
that  an  emission  occured  somewhere  in  the  tube  identified  by  the  dotted  lines.  The 
scanner  data  thus  consists  of  the  total  number  of  annihilation  counts  associated  with 
every  such  tube. 

The  discussion  can  be  made  more  precise  by  introducing  some  notation.  As- 
sume that  a unit  radius  scanner  ring  encloses  the  region  of  interest.  Let  X{r,cj)),  r G 
[0, 1],0  G [0,27t)  be  the  compound  density  function  in  polar  coordinates.  (Note  that 
A(r,  0)  = 0 for  points  between  the  patient  and  the  scanner  ring.)  Let  D be  the  (even) 
number  of  detectors.  Fix  a detector  and  call  it  detector  do.  Label  the  remaining 
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detector  ring 


detector 


detectors  di, . . . in  counter-clockwise  order.  A tube  is  identified  by  a pair  of 

detectors.  Let  T be  the  set  of  all  tubes,  T = {{di,dj}  : i ^ j,0  < i,j  < D}.  In  this 
notation,  the  data  is  a vector  v indexed  over  T ; the  element  Vt  refers  to  the  number 
of  detections  in  tube  t £ T.  The  mathematical  challenge  is  to  use  v to  estimate  the 
absorption  density  function 

Unfortunately,  the  vector  of  detections  can  be  corrupted  by  several  sources  of 
error.  These  errors  can  be  either  false  detection  of  an  emission  or  failure  to  detect 
an  actual  emission.  A false  detection  arises  when  two  photons  arising  from  distinct 
annihilation  events  strike  detectors  nearly  simultaneously.  Another  type  of  false 
detection  can  occur  when  the  line  of  flight  of  a photon  is  randomly  changed  due  to 
Compton  scattering.  An  emission  will  fail  to  be  detected  if  it  is  absorbed  by  the 
body.  Methods  to  partially  correct  for  this  effect  are  discussed  by  Ying-Lie  [41]  and 
Censor  et  al.  [7]. 
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Further  information  regarding  the  4 stages  of  PET  is  contained  in  the  survey 
articles  by  Anderson  [1],  Ter-Pogossian  et  al.  [35],  Brownell  et  al.  [5]  and  Fessler  and 
Ollinger  [29] . 

1.2  The  Standard  PET  model 

The  Shepp-Vardi  [33,  38]  mathematical  model  for  PET  models  the  tube  counts, 
Vt  as  Poisson  variables  with  mean 

n 1 

Vt=  / Q{r,4>-,t)X{r,(f))r  dr  d(f),  (El) 

Jo  Jo 

where 


Q{r,  (j)\  t)  = probability  that  a line  through  the  point  with  polar 
coordinates  (r,  4>)  and  a uniformly  random 
orientation  intersects  tube  t. 

As  stated  earlier,  the  mathematical  problem  is  to  use  the  data  v to  estimate  A(r,  (p), 

given  the  relationship  described  by  equation  1.1. 

Note  that  the  definition  of  Q{r,  0;  t)  implicitly  assumes  that  the  photons  shoot 

off  in  opposite  directions.  We  know  from  the  above  discussion  that  this  is  a slight 

simplification  of  reality.  It  is  also  apparent  that  the  Shepp-Vardi  model  does  not 

distinguish  between  the  positron  emission  process  and  the  annihilation  process.  Once 

again  we  know  from  the  earlier  discussion  that  emitted  positrons  actually  travel  a 

short  distance  before  they  collide  with  an  electron.  Thus,  the  assumption  that  the 

emission  and  annihilation  processes  are  the  same  limits  the  accuracy  of  the  model. 
1.2.1  Eourier  Solution  Methods 

Fourier  methods  rely  on  several  approximations: 

• The  detectors  are  infinitely  narrow  (i.e.  each  detector  is  a point;  each  tube  is  a 
line), 
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• all  emissions  within  a tube  have  the  same  probability  of  detection  (i.e.  Q{r,  (t>\  t) 
is  constant  for  (r  cos  0,  r sin  0)  in  tube  t), 

• the  data  is  Gaussian  (for  large  means,  Poisson  Gaussian). 

Under  these  assumptions,  the  model  reduces  to  transmission  tomography.  Hence, 
algorithms  to  recover  A(r,  0)  can  be  based  on  either  the  formula  for  inverting  the 
Radon  transform  or  the  Projection-Slice  Theorem.  These  algorithms  are  well  known 
and  described  in  Deans  [9],  Natterer  [28]  and  Lewitt  [23]. 

The  underlying  assumptions  suggest  that  the  performance  of  Fourier  methods 
will  degrade  as  the  actual  detector  size  grows.  Such  behavior  is  significant  since 
detector  size  is  limited  by  technology  and  also  because  larger  detectors  should  register 
more  counts  and  hence  provide  better  estimates  of  the  mean  number  of  emissions 
emanating  from  a tube.  Furthermore,  Fourier  methods  perform  poorly  when  the 
number  of  detections  is  low  (see  Anderson  [1]).  This  low-count  case  is  of  particular 
clinical  interest  since  it  arises  from  the  desired  scenario  of  low  radiation  doses  and 
short  scan  times.  Nevertheless,  Fourier  algorithms  are  used  in  practice  because  they 
are  fast. 

1.2.2  Solution  Methods  based  on  Discretization 

Shepp  and  Vardi  [33,  38]  presented  a statistical  approach  for  solving  model  1.1 
that  does  not  incorporate  the  simplifications  needed  for  Fourier  solution  methods. 
Their  approach  to  estimating  A(r,  0)  using  equation  1.1  begins  with  the  assumption 
that  A(r,  4>)  is  a step  function  with  respect  to  some  grid.  Let  G be  the  number  of 
grid  boxes.  The  indicator  function  for  box  i is 


1 (r  cos  (f),  r sin  0)  G grid  box  i 
0 otherwise 


The  density  function  can  be  rewritten  using  indicator  functions  as 


G 


4>)  = J^AiL(r,0). 


i=0 


(1.2) 
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In  this  case  the  model  equation  1.1  reduces  to 

j / Q{r,(f);t)li{r,(l))rdrd(f).  (1.3) 

i=i  'fo  Jo 

Thus,  i)  is  a vector  of  Poisson  variables  whose  mean  is  given  by  the  finite  dimensional 
matrix  equation 

QX  = v,  (1.4) 

where  the  elements  of  Q are  given  by 

n27T  pi 

Qti=  / Q{r,(j);t)li{r,(j))rdrd(j).  (1.5) 

Jo  Jo 

A variety  of  iterative  methods  can  be  based  on  equation  1.4.,  which,  as  the  grid 
size  approaches  zero,  is  an  accurate  physical  model  that  incorporates  finite  detector 
width.  They  perform  well  when  the  number  of  detections  is  low  (see  Anderson  [1]), 
and  are  expected  to  outperform  Fourier  methods  when  applied  to  configurations 
with  wide  detectors.  The  primary  drawback  associated  with  such  methods  is  their 
computational  complexity,  which  is  a consequence  of  the  prohibitive  size  oi  Q e 
^ standard  reconstruction  problem  could  be  based  on  G = 128^  and  D = 128 
(recall  that  D is  the  number  of  detectors),  which  would  lead  to  a matrix  that  has 
128^  elements.  Even  with  sparse  matrix  techniques  and  symmetry  considerations, 
methods  based  on  equation  1.4  are  slow. 


ML-EM  Vardi  et  al.  [38]  proposed  that  the  density  vector  in  equation  1.4 
be  determined  by  maximum  likelihood  estimation.  Since  the  data  are  independent 
Poisson  variables,  the  likelihood  function  is 


L(A)  = P(«|A)  = ne 


-vt 

Vt  ' 


(1.6) 


Dempster  et  al.’s  [10]  expectation-maximization  algorithm  is  applied  to  derive  the 
following  iterative  procedure  to  maximize  T(A), 

^tQt 


\n+l  \n  \ '' 


(1.7) 
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The  iterates  converge  monotonically  to  a global  maximum.  However,  the  limit  de- 
pends on  the  starting  guess 

Other  Iterative  Methods  Although  the  maximum  likelihood  approach  is  the 
standard  method  for  estimating  A in  equation  1.4,  many  other  objective  functions 
have  been  used.  DePierro  [8]  and  Geman  and  McClure  [13]  describe  maximum  a 
posteriori  (MAP)  methods  which  maximize  L(A)  -|-  /(A),  where  /(A)  is  the  density 
distribution  function.  The  estimate  of  / and  the  added  computational  complexity 
of  the  iteration  are  MAP’s  primary  drawbacks.  A variety  of  least  squares  (LS)  and 
weighted  least  squares  (WLS)  methods  are  discussed  in  Kaufman  [21],  Mair  et  al.  [26] 
and  Fessler  [11].  Byrne  [6]  introduced  a penalized  maximum  entropy  method  for 
which  the  algebraic  reconstruction  technique  (see  Gordon  et  al.  [15])  is  a special  case. 
It  should  be  emphasized  that  all  of  these  methods  are  iterative  schemes  based  on  the 
discrete  model. 

1.3  A New  Approach  to  Modeling  PET 

The  goal  of  this  work  is  to  develop  a fast  reconstruction  algorithm  that  does 
not  require  simplifying  assumptions  of  the  PET  model.  In  other  words,  we  seek  an 
algorithm  that  has  the  philosophical  advantage  of  incorporating  finite  detector  width, 
but  without  the  practical  speed  disadvantage  that  current  algorithms  exhibit.  To 
this  end,  we  describe  a transformation  of  the  PET  data  which  permits  a convenient 
mathematical  model  based  on  the  same  assumptions  as  the  original  Shepp-Vardi 
model.  The  following  chapters  will  analyze  this  model  and  use  it  to  develop  a new 
method  of  reconstruction. 

Transforming  the  data  Every  group  of  contiguous  detectors  identifies  an  arc. 
Let  A be  the  set  of  all  such  arcs, 

^ {{^iy  mod  £>)•••)  mod  D }i  0 <i  < D,0  < k < D}. 


(1.8) 
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Define  a 0 — 1 matrix  M indexed  over  A x T by 

Ma,t  = 


1 ant / 0 

0 otherwise 


(1.9) 


Simple  combinatorial  arguments  show  that  a scanner  ring  with  D (assume  that  D is 
even)  detectors  will  provide  data  for  ffhilhD  distinct  tubes.  The  same  scanner 
ring  will  define  {D  — 1)D  + 1 arcs.  Thus,  M is  not  invertible  since  it  is  not  even 
square.  It  does,  however,  have  full  rank. 


Proposition  1 The  transformation  matrix  M defined  by  equation  1.9  is  of  full  rank. 

Proof:  Choose  v such  that  Mv  = 0.  It  will  be  shown  that  u = 0.  The  general 
approach  is  to  show  that  Vt  = 0 for  tubes  of  the  form  t = {di,di}.  It  will  then  be 
shown  that  Vt  = 0 for  tubes  of  the  form  t = {di,di+i}.  The  process  can  be  repeated 
for  tubes  of  the  form  t = {dj,di+2}  and  so  on.  To  begin,  it  is  evident  that  the  arc 
a = {do, . . . ,c?o_i}  satisfies  aCit^tliMtET  because  it  is  the  entire  scanner  ring. 
Hence,  Ma,.  is  a row  vector  of  ones  and  so 

t t 

Now  consider  the  arc  a = {do,  • ■ • , dj_i,  . . . , do-i)  that  consists  o^  D - I detec- 
tors. It  is  clear  that  a n t 7^  0 for  every  tube  t except  (dj,  d^}.  It  follows  that 

t:^{di,di}  t 

Combine  equations  1.10  and  1.11  to  see  that 

V{d„di}  = X ^*  = 0-  (1-12) 

t tji{di,di} 

This  shows  that  uj  = 0 when  t = {di,dj,  0 < z < T>.  Next,  consider  the  arc 
that  consists  of  D - 2 detectors,  a = {do, . . . , di_i,  dj+2,  • • • , do-i}.  It  is  clear  that 
ant  7^  0 except  for  tubes  t in  the  set  S - {{d^,  d^},  {d^+i  d,  d^+i  ^od  d},  {di,  d,+i}}. 
It  follows  that 

X X 

t^s  t 


(1.13) 
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Combine  equations  1.10,  1.12  and  1.13  to  see  that 

^{di,di+i}  ^ ^di,di  ^di  + l mod  D,rfi+1  mod  D 

t t^S 

This  shows  that  = 0 when  t = {di,  dj+i  mod  £>},  0 < i < D.  Repeat  this  process 
D — 2 more  times  to  show  that  v is  identically  zero.  ■ 

We  will  model  the  transformed  data 


b = Mv.  (1.15) 

This  is  a departure  from  the  original  Shepp-Vardi  approach,  which  directly  modeled 
the  untransformed  data. 


The  mathematical  model  Note  that  ba  = {Mv)a  = Yl{teT-tna^(li}  is  the 
number  of  emissions  with  at  least  one  photon  detected  in  arc  a.  Therefore,  ba  is  a 
random  variable  with  mean 


where 


ba  = 


P{r,  (j)]  a)X{r,  4>)r  dr  dcf), 


(1.16) 


P(r,  (f)]  a)  = probability  that  a line  through  the  point  with  polar 
coordinates  (r,  (/>)  and  a uniformly  random 
orientation  intersects  arc  a. 

Since  the  Shepp-Vardi  model  is  physically  accurate,  the  purpose  of  modeling  arc  data 
instead  of  tube  data  is  not  to  introduce  a different  model.  Rather,  the  purpose  is  to 
introduce  an  equivalent  model  for  which  the  kernel  takes  a particularly  convenient 
form.  Chapter  2 will  make  this  statement  more  precise. 

Proposition  2 The  Shepp-Vardi  model  and  the  model  expressed  by  equation  1.16 
are  equivalent  in  the  sense  that 
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1.  mean  tube  detections  uniquely  define  mean  arc  detections, 

2.  mean  arc  detections  uniquely  define  mean  tube  detections. 


Proof:  These  assertions  rest  on  the  observation  '^^Ma,tQ{r,(f)-,t)  = P{r,(j)-,a).  This 
implies  that 


^27T 

= / / P{r,  (j)]  a)X{r,  (f))r  dr  d(f) 

Jo  Jo 

= / / ^Ma,tQ{rA',t)X{r,(j))rdrd(j) 

Jo  Jo  ^ 

p2n  pi 

= / Q{r,fi\t)X{r,(f))rdrd(j) 

t Jo  Jo 


ba 


= Ma,tVt- 

t 

Thus,  given  a vector  v of  tube  means  that  is  consistent  with  the  Shepp-Vardi  model, 
the  vector  of  arc  means  Mv  is  consistent  with  equation  1.16.  Conversely,  given  a 
vector  of  arc  means  b arising  from  equation  1.16,  there  exists  a vector  of  tube  means 
that  is  consistent  with  the  Shepp-Vardi  model.  Uniqueness  follows  from  the  fact  that 
M is  of  full  rank  (Proposition  1).  ■ 

Of  course,  there  are  infinitely  many  arcs  other  than  the  ones  for  which  we  have 
data.  The  set  of  all  arcs  can  be  parameterized  by  {p,9),—l  < p < 1, 0 < 9 < 2tt.  We 
identify  the  pair  (p,  9)  with  the  arc  Op^e  defined  by  moving  counter-clockwise  from 
endpoint  ci  to  endpoint  C2,  where 


ei{p,9)  = (pcos0 -h  \/w^p2sin0,  psin0  - -v/W-p^cos^),  (1-17) 

^2(9, 9)  = {pcos9  - \/l  - p2sin0,psin0  + y/l-  p‘^cos9).  (1-18) 

This  parameterization  is  illustrated  in  Figure  1.2.  The  arc  complementary  to  the 
bold  arc  identified  in  the  figure  would  have  parameters  {9  + tt,  -p).  Every  point  on 
the  unit  disk  is  similarly  associated  with  two  arcs.  This  explains  the  range  of  the 
parameters  p and  9. 
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arc 


Even  though  we  only  have  data  for  finitely  many  arcs,  we  consider  a general 
model  in  which  there  is  a continuum  of  arcs 


where 


^27T  pi 

Kp,^)=  / P{r,(f);p,e)X{r,(j))r  drdcf), 

Jo  Jo 


(1.19) 


b{p,  9) 

P{r,  0;  p,  6) 


= mean  number  of  emissions  with  at  least  one  photon 
detected  in  arc 

= probability  that  a line  through  the  point  with  polar 
coordinates  (r,  0)  and  a uniformly  random 


p,e- 


orientation  intersects  arc  a 


CHAPTER  2 

ANALYSIS  OF  THE  MODEL 


This  chapter  examines  the  arc-based  model  for  PET  introduced  in  Chapter  1. 
Section  2.1  begins  the  discussion  by  introducing  orthogonal  series  expansions  for 
and  A(r,  0).  The  particular  orthogonal  functions  that  are  used 
are  Bessel  functions,  trigonometric  functions  and  Tchebyshev  polynomials.  When 
these  series  expansions  are  substituted  into  the  model,  it  reduces  to  a set  of  infinite- 
dimensional matrix  equations.  These  matrix  equations  are  recognized  as  a potential 
basis  for  a new  class  of  PET  reconstruction  algorithms.  Such  an  approach  requires  a 
framework  which  facilitates  study  of  the  matrix  structure  and  computation  of  the  ma- 
trix elements.  This  motivates  the  deeper  analysis  which  is  the  subject  of  section  2.2. 
In  that  section,  it  is  discovered,  that  the  PET  kernel  is  superharmonic  and  its  Riesz 
Decomposition  is  derived.  The  remainder  of  the  chapter  uses  this  decomposition  to 
find  explicit  formulae  for  the  series  coefficients  of  P(r,  (/»;  p;  9). 

2.1  Model  Formulation  in  terms  of  Orthogonal  Functions 

This  section  will  relate  series  expansions  for  \{r,4>),P{r,(l);  p,9)  and  b{p,6). 
Before  this  problem  is  tackled,  a brief  review  of  the  orthogonality  relations  for  Tcheby- 
chev  polynomials  and  Bessel  functions  may  prove  useful. 

Tchebyshev  Polynomials  The  Tchebyshev  polynomial  of  degree  m is  con- 
ventionally denoted  by  Tm{p)  and  is  defined  by 

T’m(p)  = cos(mcos“^(p)),  -l<p<l.  (2.1) 


12 


13 


Trigonometric  identities  can  be  applied  to  show  that  T^(p)  is  in  fact  a polynomial  of 
degree  m that  satisfies  the  recurrence 

T,{p)  = 1 
Trip)  = p 

Tm+\{p)  = “^pTmip)  - Tm-l{p),  m > 0. 


(2.2) 


They  also  satisfy  the  orthogonality  relation 

Up)T,{p) 


L 


dp 


0 z 7^  j 

f * = jV  0 , 

7T  Z = j = 0 


(2.3) 


•yi  - 

The  following  Theorem  regarding  series  expansions  follows  directly  from  the  well 
known  analogue  for  trigonometric  series  (See  Littlewood  [24]). 

Theorem  3 (Tchebyshev  Series  Expansion)  ///  G £^[0,27r]  then 

OO 

/(P)  = ^fmTm{p), 

where 

fm  - \ y rl  r,  X Tm{p) 


m=0 


Proof:  Define  a function  g G £^[0,  27r]  by  g{6)  = /(cos(6>)).  It  is  well  known  that 
the  trigonometric  series 

OO 

fmCos{me), 

m=0 

where 

f = / A -jf  ^ cos(mP)  dp,  m = 0 

^ \ ^ /o  ^ otherwise  ’ 

converges  to  g in  £^[0,27t].  The  theorem  follows  from  the  change  of  variable  6 = 

cos~^  p.  ■ 
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Bessel  Functions  The  Bessel  function  of  the  first  kind  of  order  > — | has 
the  series  representation 


= ’,y 

2 ^j!r(A:  + j + l) 


(2.4) 


Like  Tchebyshev  polynomials,  Bessel  functions  also  satisfy  a recurrence  formula  and 
an  orthogonality  relation.  The  recurrence  is 


2kJk(r)  = rJk-i(r) +rJk+i(r),  k > 1. 


(2.5) 


The  orthogonality,  however,  is  over  the  interval  [0, 1]  with  respect  to  the  weight  r, 
'1 


n ^ m 


I Zji  k)  dv  — ^ 1 t2  / \ ) 

Jo  I ^Jk+ii^n,k)  n = m ' 


(2.6) 


where  Zn,k  is  the  positive  real  root  of  the  Bessel  function  of  order  k,  Jk{r).  The 
following  four  indefinite  integrals  involving  Bessel  functions  will  be  useful  later  in  this 
discussion  and  can  be  found  in  Bowman  [4].  For  any  real  z and  k > 0, 

f k+i  T /■  \ t ^^^^Jk+iirz) 

/ r'^^^Jk{rz)dr  = 

J r~'^+^Jk{rz)dr  = 


+ C 


fc-i 


J r\og{r)Ja{r)dr  = Jo(r)  + ’'log(r)Ji(r)  + C 
J r log(r)  dr  = + r log{r)  J,{r2)  + C. 


(2.7) 

(2.8) 
(2.9) 

(2,10) 


Note  that  equation  2.10  is  a direct  consequence  of  equation  2.9.  Theorem  5 identifies 
functions  that  can  be  expanded  in  Bessel  series.  It  makes  use  of  the  following  result 
from  Watson  [40]. 


Theorem  4 (Uniform  Convergence  of  Bessel  Series)  Let  \/{r)f{r)  G iC^O,  1]. 
In  addition,  let  f be  continuous  on  [0, 1]  and  of  bounded  variation  in  (5, 1 — 4).  The 
Bessel  series  expansion  of  f{r)  converges  to  f{r)  on  [0, 1]  and  the  convergence  is 
uniform  on  (4, 1 — 4). 
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Theorem  5 (Bessel  Series  Expansion)  Let  k > \.  If  f e £^[0, 1]  then 

00 

/(^)  = ^fnJk{Zn,kr), 

where 


n=l 


fn  = 


Jl+M 


[ f{r)Jk{zn,kr). 

n,k)  Jo 


Proof:  We  need  to  show  that  the  Bessel  functions  are  complete  in  £^[0, 1].  Since  the 
continuous  functions  of  bounded  variation  in  £^[0, 1]  are  dense  in  £^[0, 1],  it  sufhces 
to  show  that  for  every  /(r)  G C[0, 1]  n £^[0, 1]  of  bounded  variation  and  for  every 
e > 0,  there  is  a finite  combination  of  Bessel  functions  s(r)  such  that  ||s  - /H2  < e. 
To  achieve  this,  reinterpret  the  orthogonality  with  respect  to  weight  r as  ordinary 
orthogonality  of  the  functions  y/rJk{rZn,h)-  Similarly,  the  Bessel  expansion  of  /(r)  has 
the  same  coefficients  as  the  expansion  of  y/rf{r)  with  respect  to  the  orthogonal  system 
VrJk{rZn,k)-  f{r)  G £^[0, 1]  implies  that  ^/ff{r)  e £^[0, 1].  Since  £^[0, 1]  C £^[0, 1], 
we  have  \/r/(r)  G £^[0, 1].  Therefore,  Theorem  4 implies  that 

00 

y/{r)f{r)  = J^c„x/(r)Jfc(r2„  jfc),  r G (0,1). 

71=1 

Bessel’s  inequality  (see  Rudin  [32])  implies  that 


00  > 


pi  pi 

/ rf{r)  dr>^cl  / rJl{rZn,k)  dr  = L. 
•^0  „=i  ^0 


It  follows  from  the  finiteness  of  this  sum  that  there  exists  an  N such  that 

N 1 

/ rJl{rZn,k)  dr\  < e. 

n=l  Jo 

Since  the  terms  in  the  summation  are  nonnegative, 

00 

XI  4 / rJl{rZn,k)  dr  < e. 
n=AT+l  Jo 

Equation  2.11  shows  that  the  partial  sums  converge  to  /(r)  since 

00  00 

ll/(^)  ~ 'y  ^ ~ II  ^ ] c„\/r  J/j(rz„  *,)  H2  < ^ / f Jk{'^^n,k)dr  < c. 

n-1  n=N+l  n=N+l  Jo 

(2.12) 


(2.11) 
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Double  and  Triple  Series  Expansions  We  conclude  this  section  by  combin- 
ing Tchebyshev,  Bessel  and  trigonometric  series  as  described  in  the  Theorems  6,7 
and  8. 


Theorem  6 (Trigonometric-Tchebyshev  Series)  If  the  function  f{p,  6)  is  in 
L^([— 1, 1]  X [0,  27t])  then 


fc=— oo  m=0 

where  the  coefficients  are  given  by 

-1  r27r 


/I  rZTT 

J 


\/T^ 


m = o ■ 


Theorem  7 (Trigonometric-Bessel  Series)  If  the  function  f(r,6)  is  in 
L^([0, 1]  X [0,  27t])  then 

00  OO 
k=—oo  n=l 

where  the  coefficients  are  given  by 


ni  pZTT 

Theorem  8 (Trigonometric-Bessel-Tchebyshev  Series)  If  the  function 
f{r,  6,  p)  is  in  1]  x [0,  27t]  x [-1, 1])  then 

OO  OO  OO 

f{r,e,p)=  E J2T, , 

k=—cx3  n=l  m—0 

where  the  coefficients  are  given  by 

'1  />!  r2-K 


/i  pi  pZTT 

I I /(^)  ^) 

1 Jo  Jo 


P) 


'^m{p)  -I 


ikO 


J\k\{rzn,\k\)r  dedpdr, 


^k,Ti 


, m 7^  0 


-ji-  , m = 0 
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Application  to  PET  Recall  that  the  PET  model  in  polar  coordinates  is 


(2.13) 


The  variables  in  equation  2.13  lie  in  the  intervals  0<r<l,-l<p<l  and 
0 < 0 < 27t.  Note  that  these  intervals  correspond  precisely  to  the  intervals  over  which 
the  Bessel  functions,  Tchebyshev  polynomials  and  trigonometric  functions  exhibit 
their  respective  orthogonality.  Furthermore,  the  weight  r that  arises  from  the  change 
to  polar  coordinates  also  happens  to  be  the  weight  function  needed  for  the  Bessel 
orthogonality  relation  described  in  equation  2.6.  These  observations  suggest  that  we 
should  examine  the  series  expansions 


where  A;  G Z,  n € N,  m G Nq.  Here  the  natural  numbers  are  denoted  by  N,  the  noneg- 
ative integers  are  denoted  by  Nq  and  Z refers  to  the  integers.  Note  that  equation  2.16 
uses  the  circular  symmetry  of  P{r,  (f);  p,  6)  that  is  a consequence  of  the  scanner  geome- 
try. It  will  be  shown  later  in  this  chapter  the  P{r,  0;  p,  6)  is  continuous  and  bounded, 
hence  the  expansion  above  is  valid.  It  follows  that  the  above  expansion  is  also  valid 
for  b{p,6)  so  long  as  A(r,  0)  G £^[0,1]  x [0,  27t].  The  PET  model  is  now  interpreted 
in  terms  of  the  coefficients  in  these  series  expansions. 

Proposition  9 The  coefficients  in  the  trigonometric- Tchebyshev  expansion  ofb{p,  9) 
are  given  by 


(2.14) 


k,n 


k,m 


k,m,n 


(2.17) 


n 


18 


Proof:  Substituting  the  series  expansions  for  A and  P into  the  PET  model  in  equa- 
tion 2.13  gives 


This  equation  can  be  rearranged  as  follows  so  that  the  role  of  orthogonality  is  evident: 


The  orthogonality  of  the  trigonometric  functions  and  that  of  the  Bessel  functions 
serve  to  simplify  this  expression  to 


Theorem  6 shows  that  this  is  the  trigonometric-Tchebyshev  series  expansion  of  6(p,  9). 

■ 

This  tells  us  that  for  every  k,  Afc_.  satisfies  an  infinite  dimensional  matrix  equa- 
tion. Chapter  3 uses  this  matrix  equation  as  the  starting  point  for  a new  class  of 
PET  reconstruction  algorithms.  At  this  stage,  the  choice  of  Tchebyshev  polynomials 
is  quite  unmotivated.  However,  subsequent  sections  will  reveal  a connection  between 
Bessel  functions,  Tchebyshev  polynomials  and  the  kernel  P{r,4>]  p,9).  The  next  sec- 
tion presents  a decomposition  for  P(r,  0;  p,  9)  that  can  be  used  to  calculate  the  matrix 
elements. 


Suppose  we  are  given  an  arc  ap^g  with  arclength  ap^g.  Let  Lp^g  denote  the  line 
segment  that  joins  the  arc  endpoints  ei  and  C2.  Lp^  can  be  parameterized  as 

Lpp{t)  = {pcos9  + tsin9,  psin9  - tcos9),  -\/l~p^<t<  y/l-p2.  (2.18) 

Denote  the  subset  of  the  unit  disk  enclosed  by  the  arc  and  the  line  by  R.  Denote 
the  remaining  subset  of  the  unit  disk  by  S.  This  labeling  scheme  is  illustrated  in 


k,n  k' ,n' ,m 


2.2  Riesz  Decomposition  of  the  PET  Kernel 
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Figure  1.2  and  will  aid  in  the  proof  of  several  propositions  that  follow.  Associated 
with  each  arc  is  an  angle-of-view  function  defined  by 

fi*)  {x,y)€S 


Qp,e{r,(f))  = 


cos" 


cos 


-1 


(2.19) 


where 


y = {r  cos  4>,  r sin  4>)  — el{p,  9) 
z = (r  cos  0,r  sin  (;;!»)  — e2(p,0). 

As  the  name  suggests,  the  angle-of-view  function  has  a simple  geometric  interpreta- 
tion as  the  angle  between  vectors  formed  by  joining  a point  in  the  unit  disk  to  the 
endpoints  of  an  arc. 

Proposition  10  P is  the  minimum  of  two  functions, 

P{r,  0;  P,  9)  = min{l,  ^Qp,e{r,  </>)}.  (2.20) 


Proof:  Simple  geometric  arguments  show  that  P(x;  p,  9)  can  be  computed  as 


<f>\  P,  0) 


1 (r  cos  (/),  r sin  0)  € R 

(r  cos  r sin  G S' 


The  proposition  follows  from  the  observation  that  0p,e(r,  0)  > tt  for  (r  cos  0,  r sin  4>)  e 
R.  m 


The  representation  in  Proposition  10  gives  immediate  insight  into  important 
properties  of  the  kernel,  as  evidenced  by  Corollaries  11  and  12. 


Corollary  11  P{r,  (j)\p,9)  is  non-negative  and  superharmonic  in  {r,(j)). 

Proof:  The  harmonic  measure  associated  with  arc  Op^g  is  the  function  that  satisfies 

Aujp^g(x)  = 0 in  the  unit  disk,  along  with  the  boundary  conditions 

. ^ _ J 1 if  X lies  on  arc  Op^g 

\ 0 if  X lies  on  the  arc  complimentary  to  Op^g 
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We  can  write  the  angle-of-view  from  a point  to  an  arc  in  terms  of  harmonic  measure 
as 


QpA'f'^  A 


TTLUp^g{r  COS  (f),  r sin  (/>)  + 


2 ■ 


It  is  therefore  evident  that  equation  2.20  expresses  P as  the  minimum  of  harmonic 
functions.  ■ 


Corollary  12  P admits  a Reisz  Decomposition  as  the  sum  of  a Green’s  potential  of 
a measure,  Pp^g,  and  a harmonic  function,  H[--,p,0), 


p{r,<t>;  P,  = Gpff(r,  ^)  + H(T,<j>\  p,  0) . 


(2.21) 


In  this  decomposition,  G is  the  Green’s  function  for  the  unit  ball,  which  is  defined  as 

G(a:;!/)  = < -log(||!,||)  z€B\{x},i  = 0,  (2.22) 

t oo  y = X 

where  x*  = is  the  inverse  of  x with  respect  to  the  unit  circle.  Gpp^g{r,(j))  satisfies 

GppAfA)  = / G{x-,  p cos  6,  p sin  9)  dpp Ax),  (2.23) 

Jb 


where  B is  the  unit  disk. 


Proof:  This  is  the  Riesz  Decomposition  Theorem  (see  [19]).  ■ 

The  remainder  of  this  section  is  devoted  to  an  examination  of  the  components 
of  the  Riesz  Decomposition  of  P(r,  </>;  p,  6).  Proposition  13  characterizes  the  harmonic 
part  in  terms  of  the  Poisson  Integral  Formula.  The  succeeding  proposition  turns  our 
attention  to  the  Green’s  potential.  Together  they  provide  expressly  stated  formulas 
for  the  Riesz  Decomposition  that  will  be  used  to  calculate  the  series  coefficients  in 
section  2.3. 

Proposition  13  P can  be  expressed  as 


(2.24) 
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where 


hp{t) 


1 — cos  ^ p <t  < cos  ^ p 

^ otherwise. 

7T 


(2.25) 


Proof:  Referring  to  the  Riesz  Decomposition  in  equation  2.21,  H is  the  function  that 
is  harmonic  on  the  unit  disk  and  shares  the  the  same  boundary  values  as  P.  The 
circular  symmetry  of  P implies  that  H{r,  (j)\  p,  9)  = H{r,  (j)  - 6]  p,0).  Any  harmonic 
function  on  the  unit  disk  can  be  expressed  in  terms  of  its  boundary  values  by  means 
of  the  Poisson  Integral  Formula  (see  Rudin  [32]).  The  Poisson  integral  representation 
for  H (in  polar  coordinates)  is 


H{r,(f);p,9)  = H(r,(f)  - 9]  p,0) 


1 r 


1 — r^ 


2r  cos(0  ~ 9 — t)  + r'^ 
where  equation  2.20  gives  the  boundary  values 


hp{t)  dt, 


(2.26) 


hp{t) 


1 if  lies  on  arc  (p,  0) 
^ otherwise 


(2.27) 


The  series  expansion  of  the  Poisson  Kernel  (see  Rudin  [32])  is 


1 — 2r  cos{4>  — 9 — t)  + r"^ 


(2.28) 


Substituting  the  series  expansion  of  the  Poisson  Kernel  into  equation  2.26  gives  the 
representation  in  the  proposition.  ■ 

As  noted  earlier,  the  next  proposition  addresses  the  Green’s  potential  of  the 
measure  in  the  Riesz  decomposition  of  P{r,  0;  p,  9). 


Proposition  14  The  Green’s  potential  of  Pp^g  in  Corollary  12  can  be  expressed  as 


Gpp,e{r,  (j)) 


2 

= — / G(r  cos  0,  r sin  (/>;  pcos0  + t sin0,  psin0  — t COS0)- 

7T  1 


p^  — p 


(St29) 
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Proof:  This  proof  abuses  notation  slightly  by  referring  to  the  PET  kernel  in  Carte- 
sian coordinates.  By  P{x-,p,9)  we  mean  P{r,(f)\  p,0),  where  x = (r  cos  (/),  r sin  (/>). 
This  should  not  cause  confusion  since  this  proof  makes  exclusive  use  of  Cartesian 
coordinates.  In  addition,  the  usage  is  clear  from  the  number  of  arguments.  The  same 
convention  is  followed  for  H{x-,p,6)  and  Ppp{x). 

Observe  that  applying  the  Laplacian  to  P gives  the  result  AP{x-,p,9)  = 
^Gppfi{x)  + AH{x-,p,9)  — -ppfi{x).  We  interpret  this  in  the  sense  of  distributions 
to  mean  that  for  every  test  function  T, 

f {AP{x-p,9))'if{x)  dA{x)  = [ P{x-,p,9)A^{x)  dA{x)  = - [ T(x)  dpp^g{x). 

Jb  Jb  Jb 

(2.30) 

We  may  rewrite  this  integral  over  the  disk  as  an  integral  over  the  disjoint  regions  R 
and  S, 

[ P{x-,p,9)A'i>{x)dA{x)=  [ P{x]p,9)A'ii{x)dA{x)+  [ P{x]  p,9)A'^{x)  dA{x) 
Jb  Jr  Js 

(2.31) 

Green’s  second  identity  (see  [19])  can  be  applied  to  each  integral  to  obtain 


IR 


P{x]  p,9)A^{x)  dA{x)  = 


f r 8P 

/ (AP(x;  p,  0))T(a;)  dA{x)  + / P(x;  p,  9)^{x)  - T(x)  — (a:;  p,  9)  da{x) 

Jr  JdR  our  driR 

f P(x;  p,  0)AT(x)  dA{x)  = 

Js 

f /■  (9P 

/ {AP{x-,  p,9))^{x)  dA{x)  + / P{x-,p,9)-—{x)-'ij{x)j—{x-,p,9)da{x) 

Js  Jds  (J'^s  ons 

Since  P is  harmonic  in  region  R and  in  region  S,  Green’s  identity  simplifies  to 

j P{x-p,9)A<^{x)dA{x)  = / P{x-p^9)  — {x)-'^{x)  — {x\p,9)da{x) 

Jr  JdR  obr  our 

f r BP 

P{x-p,9)A'^{x)dA{x)  = / P{x-p,9)^{x)-<l>{x)^{x;p,9)da{x) 

Js  Jas  (JBs  Bus 
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9)  da{x) 


Furthermore,  since  ^ is  a test  function  whose  support  is  a compact  subset  of  the  unit 
disk,  = 0 when  |x|  = 1.  Thus,  we  can  change  the  region  of  integration  to  obtain 

/■  r 8P 

P{x-p,9)A'^{x)dA{x)  = / P{x-,p,9)- — (x)  - v[r(x)- — {x]  p,9)  da{x) 

Jr  Jbp^g  oRr 

f f BP 

P{x;p,9)A'^{x)  dA{x)  = / P{x- p,9)-—{x)  - ^{x)^{x]  p 

Js  Bus  Bus 

We  may  now  substitute  this  application  of  Green’s  second  identity  into  equation  2.31. 
Noting  that  ^{x)  = --^{x)  when  x is  on  line  segment  Lp^e,  we  have  the  following 
equation, 

J^P{x-p,9)A^{x)dA{x)  = -J^  (^^{x)^{x-,p,9)  + ^{x)^{x- p,9)^  da{x). 

(2.32) 

Since  P is  constant  on  R,  the  last  integrand  in  equation  2.32  is  zero.  Therefore,  we 
may  rewrite  equation  2.32  in  terms  of  harmonic  measure  as 

[ P{x-,p,6)A^{x)dA{x)  = -[  'i>{x)^~^{x)  da{x).  (2.33) 

Jb  Jl^^s  Bus 

Referring  to  the  parameterization  of  the  line  Lp^,  the  normal  derivative  of  harmonic 
measure  can  be  computed  using  trigonometric  arguments  to  obtain 

{p  cos  9 + tsin9,  psin9  — tcos9)  = — n ^ n- 
Bus  ’ -K  I-  p^  -P 

By  equations  2.30  and  2.33,  this  leads  to  the  following  characterization  of  the  measure 
Pp,e'- 


f 2 f\  2 

/ 'i>{x)  dpp^g{x)  = - 'i’{pcos9  + tsin9,  psm9  - tcos9)— — dt. 

Jb  1-  p^  -P 


As  a consequence  of  the  above  propositions,  the  Riesz  Decomposition  is  now 
clearly  defined.  The  corollary  below  simply  summarizes  the  analysis  in  this  chapter; 
the  PET  model  is  rewritten  in  terms  of  the  Riesz  Decomposition. 
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Corollary  15  The  PET  model  expressed  in  polar  coordinates  can  be  written 

Kp,^)=  / H{r,(j);p,9)X{r,(f))r  drdcj) 

Jo  Jo 


'0  Jo 

r2n  rl 


2 /•■in  pl 

H — / / / G{r  cos(f),rsm(j)\  pcosO  -\-  tsmO,  psinO  — tcosB) 

Jo  Jo  4-,/i-p2 


V 

dt  A(r,  4>)r  dr  d(j) 


1 — fp  — t'^ 

2.3  Calculations  for  the  Series  Coefficients  of  H 


This  section  and  the  following  section  are  devoted  to  finding  the  trigonometric- 
Bessel-Tchebyshev  series  coefficients  for  P(r,  p,  4>)  by  finding  the  coefficients  for  the 
potential  and  the  harmonic  function  that  arise  from  its  Riesz  Decomposition.  We 
begin  by  finding  a series  expansion  for  the  harmonic  function  H (r,  0;  p,  9)  of  the  form 

H{r,  (j)]  p,  9)  = H{r,  (p-9;p,0)=  ^ Hk^rn,nJ\k\{rZn,\k\)Tm{p)E^^‘^~^\  (2.34) 

n,kym 

where  A:  G Z,  n e N,  m e Nq.  Just  as  in  equation  2.16,  equation  2.34  uses  the 
circular  symmetry  inherent  in  H{r,(f)-,  p,9).  As  noted  in  the  proof  of  Proposition  13, 
this  circular  symmetry  arises  from  the  fact  that  H{r,  (j)]  p,9)  is  determined  by  the 
boundary  values  of  P{r.,(f)-,  p,4>),  and  we  observed  earlier  that  P{r^(f)\  p,<f)  exhibits 
circular  symmetry. 


Proposition  16  The  trigonometric-Bessel-Tchebyshev  series  coefficients  for  H are 
given  by 


Hk,m,n  — ^ 


3^1 

— 8 

n m J\k\-\-\{,/^n,\k\)^ny\k\ 
2 


,1*1 


A;  = 0,  m = 0 
A;  = 0,  m / 0 
^ 0 

A;  0,  m = 0, 
otherwise 


(2.35) 


/r'^{k'^-m'^)Jlk\+i{z„^\k\)zn,\k\ 

Proof:  We  can  apply  Theorem  8 to  the  Poisson  series  representation  of  H in 
equation  2.24  to  write  the  series  coefficients  Hk,m,n  as 

T-  L C I'  L dtdrd^dp. 

(2.36) 
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Grouping  the  terms  yields 

(2.37) 

This  expression  can  be  simplified  using  the  orthogonality  of  the  trigonometric  func- 
tions to  obtain 


rl^l+V|fc|(rz„,|,|)  dr  ['  f dtdp.  (2.38) 

Jo  J-l  J-n  y/1  — 

The  integral  involving  the  Bessel  term  can  be  evaluated  with  the  aid  of  equation  2.7 

to  obtain 

P r (()  it  10^  dp,  (2.39) 

J-,  J-, 

It  follows  from  the  definition  of  hp{t)  in  equation  2.25  that  the  innermost  integral  in 
equation  2.39  can  be  expanded  as 


f e ^^^hp{t)  dt 

J —7T 

/cos-i(p) 

' 

• COS“l(p) 


7T 


£ 


— cos  ^ (p) 


dt  + I dt  I . (2.40) 

J cos~qp) 

Thus,  we  can  apply  the  change  of  variable  x = cos~^(p)  to  equation  2.39  to  arrive  at 
('k,m,nJ\k\  + l{^n,\k\) 


^n,\k\ 


f cos{mx)  ( [ e dt+-  [ e dt  + - f e dt\  dx. 
0 —X  J — n ^ Jx  J 

(2.41) 


We  first  turn  our  attention  to  the  case  where  k ^ 0.  Computing  the  integrations 
with  respect  to  t yields 

'^Pk,m,nJ\k\+l{^n,\k\)  f 


kz. 


'n,|fc| 


/ cos(ma;)(f 

Jo 


—ikx  ^ikx  _L  Jp^ikx  ^ ^ik-K  , ^ „—ikir 


e -\ — e 

7T 


-e  H — e~ 

7T  7T 


7T 


(2.42) 


This  can  be  rewritten  using  the  identity  2i  sin(z)  = to  obtain 

‘J‘Pk,m,nJ\k\  + \{^n,\k\)  f 


kz. 


n,\k\ 


r X 

/ cos(mx)(l )sin(A:a;)  dx. 

Jo  7T 


(2.43) 
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The  following  two  facts  aid  in  the  computation  of  the  integral  in  2.43: 

/ cos{mx)  sm{kx)  dx  = 

Jo 


f 


X 

— cos{mx)  sin{kx)  dx 

7T 


0 

= k^ 

(2.44) 

- otherwise 

0 

o 

II 

II 

1 

4k 

fc(_l)*:+m 

m?—k'^ 

m2  = it  V 0 
otherwise 

(2.45) 

The  integral  in  equation  2.41  can  be  obtained  in  a straight-forward  manner  for  the 
case  k — 0.  This  gives  a convenient  formula  for  the  series  coefficients  of  H, 


CQ.n.O'Tl  (^n)  47T^ 


H; 


k,m,n 


^0,n,m  1 (^n  ) 4 


= 


A:  = 0,  m = 0 
A:  = 0,  m 0 


= k"^  0 


^k,m,n'^\k\  + li^n,\k\) 

Otherwise 

(A:  m 


(2.46) 


Substituting  the  definition  of  Ck,m,n  from  equation  8 into  equation  2.46  yields  the 
proposition.  ■ 


The  proof  of  proposition  16  reveals  a general  relationship  between  the  Poisson 
Integral  Formula  and  Bessel  functions.  The  proof  also  identifies  a special  relationship 
between  the  boundary  values  of  H{r,  (j),  p,  9)  and  the  Tchebyshev  polynomials.  Thus, 
the  orthogonal  functions  in  section  2.1,  which  may  have  appeared  rather  arbitrary, 
were  carefully  chosen  to  exploit  their  association  with  the  PET  kernel. 

A useful  restatement  of  formula  2.35  is  presented  in  Corollary  17.  We  will 
see  in  Chapter  3 that  this  restatement  has  important  computational  implications  for 
both  computing  the  coefficients  and  solving  the  system  of  equations. 


Corollary  17  77*, is  the  rank-1  matrix  uv^  where 

1 


Ur. 


( 4 

3 


= < 


— ,n  e N, 

(2.47) 

,,|fc| 

A;  = 0,  m = 0 

A:  = 0,  m 0 

m2  = A;2  / 0 , m e No. 

(2.48) 

A;  ^ 0,  m = 0, 
otherwise 
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2.4  Calculations  for  the  Series  Coefficients  of  G/, 

Next,  we  find  find  a formula  for  the  trigonometric-Bessel-Tchebychev  series 
expansion  coefficients  of  the  Green’s  Potential,  G/i.  These  coefficients  are  determined 
by 

Gfx{r,  (f)-  p,  6)  = Gp{r,  (j)-e-p,0)  = ^ Gpk,m,nJ\k\{rZn,\k\)Tm{p)e^^^'^~^\  (2.49) 

n,k,m 

where  A:  € Z,  n G N,  m € Nq.  Let  us  first  examine  the  functions 

r2-K 

Tk,m{s)=  / Tm{scos(l))e~'^'^  d(j).  (2.50) 

Jo 

The  following  three  simple  facts  regarding  these  functions  will  be  useful. 


Proposition  18 


Ti 


k,m 


^ 0 |A;|  > m 
0 k even,  m odd 


Otherwise,  T*, ; 


0 k odd,  m even 
is  a polynomial  containing  terms  of  degree  m,  m 


(2.51) 


2,  m — 4, . . . , A:. 


Proof:  Suppose  that  m is  even.  The  Tchebyshev  polynomial  Tm{x)  consists  of  terms 
of  degree  0,  2, . . . , m.  Write  T^(x)  = aQ  + a2X^  + ■ • • + The  definition  of  Tk^rn{s) 


yields 


Tk,m{s)  = / Tm{s  COS  (j))  COS  k(j)  d(j) 

Jo 

m 

r2ir  2 

U2jS^^  cos'^^  (j)  COS  k(j)  d(j) 


° j=0 


^02j  j cos^-^  (/)  cos  k(j)d<^  s^^  (2.52) 


This  shows  that  Tk^rn{s)  is  a polynomial  consisting  of  even  powers  of  s up  to  m.  We 
will  now  perform  the  integration  to  show  that  many  of  these  coefficients  are  zero. 
Repeated  application  of  the  trigonometric  identity  cos‘^{(f))  = 4(1  + cos  2^)  shows 
that  there  exist  coefficients  b,  such  that 

3 

cos^  j4>  = cos(2z0). 

i—\ 


(2.53) 
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This  implies  that 


It  is  clear  that  equation  2.54  reduces  to  zero  when  2j  < A:  or  if  A:  is  odd.  Extend 


this  result  by  substituting  equation  2.54  into  equation  2.52  to  see  that  Tk,m{S)  = 0 
if  m < A:.  This  proves  the  proposition  for  m even.  A similar  argument  proves  the 
theorem  for  the  case  where  m is  odd.  ■ 

Proposition  19  1 and  —1  are  roots  of  the  polynomial  Tfc_^(s),  m ^ k. 

Proof:  Equation  2.50  shows  that  ±1  are  roots  of  k since 


0,  m ^ k 
2tt,  m = k = 0 

7T,  m = \k\  ^ 0 


Proposition  20  The  functions  Ti-^rn{s)  satisfy  the  recurrence 


(2.55) 

(2.56) 


7TS‘ 


,k 


Proof:  We  will  make  use  of  the  trigonometric  identity 


cos(0)  cos{k(l))  = cos((A:  ± !)(/>)  ± sin((/>)  sin{k4>). 


(2.58) 


Apply  this  identity  twice  to  obtain 


2cos(0)  cos{k(j))  = cos((A;  + l)^)  + cos((A;  - 1)0). 


(2.59) 
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Substitute  the  Tchebyshev  recurrence  formula  2.2  into  the  definition  of  Tk,m  given  in 
equation  2.50  to  obtain 


Tk,m{s)  = / 2s  COS  COS  (j))  COs{k(j))  — Tm-2{s  COS  d(f) 

Jo 

p27T 

= S 2cos(f)Tm-i{scos(j))cos{k(f))d(l)-Ti,,m-2 
Jo 

Now  substitute  equation  2.59  into  equation  2.60  to  see  that 


(2.60) 


Tk,m{s) 

p2TT 

= S Tm-i{scos(j)){cos{{k  + l)(f))  + cos{{k-l)(t)))d(l)-Tk,m-2 

Jo 

= +Tk+l,m-l{s))  —Tk,m-2{s)-  (2-61) 

The  initial  values  follow  from  straight-forward  integration.  ■ 

We  are  now  in  a position  to  develop  a formula  for  the  coefficients  of  the  Green’s 
potential  of  the  measure  identified  in  the  Riesz  Decomposition.  We  begin  by  finding 
the  coefficients  of  the  Green’s  function  and  then  use  that  result  to  find  the  coefficients 
of  the  Green’s  potential. 


Proposition  21  Let  Gk,n  be  the  coefficients  in  the  trigonometric-Bessel  expansion 
of  the  Green’s  function  for  the  unit  disk, 

"1  />27T 


Gk,n{^)  = 77- T ! f G {r  COS  if,  r sin  i)- s,0)e  ^^’f’rJik\(rZn,\k\)  dip  dr. 

2t^ J\k\\Zn,\k\)  Jo  Jo 


(2.62) 


The  integrals  can  be  evaluated  to  obtain 


Gk,n^^') 


' J\k\  (^■2'7i,|fc|)  ■ 


Proof:  The  assertion  follows  from  the  power  series  expansion 


log(l  -z)  = -z-  - ^z^  - ^z*  - O(z^), 


(2.63) 


(2.64) 
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valid  for  \z\  < 1.  Suppose  0 < r < s.  From  the  definition  of  the  Green’s  function  in 
equation  2.22,  we  can  write 


G(r  cos  tp,  f sin  ip]  s,  0) 
= log 


s|i  - re^^\ 


s — re 


iip  I 


d 


log(s)  + log  I - — — log  Is  — 

s 

= log(s)  + Re  log(-  — — Re  log(s  — re®®^) 

= — log(s)  + Relog(l  - sre®®^)  — Relog(l  — -e®^). 

s 

Substitute  the  expansion  from  equation  2.64  into  equation  2.65  to  obtain 


(2.65) 


G{r  cos  xp,  r sin  ip-,  s,  0) 


= -log(s)  + -Re  J]-(rs)'®e''®®^  + ReVi  (-) 

. rC  K \ S ' 

k>0  k>0 

= -log(s)  + Re^  ^ (1  - s^*’)e®*’^ 

. K \ S ' 
fc>0 

= ~ log(^)  + ^ cos{kxp) 

k>0  ^ 

= -logM  + i$:t(d‘(l-s“)e« 

Jt^O 


1 /r\k 


(2.66) 


Observe  that  equation  2.66  is  a trigonometric  series.  Symmetry  of  the  Green’s  func- 
tion can  be  used  to  obtain  analogous  results  for  s < r.  It  follows  that 

— log(r)  A;  = 0,  s < r 

^(i)|fc|(l  _ r2|fc|)  k^0,s<r 


1 

— / G{r  cos  xp,r  sin 'ip-,s,0)e~ dip 
27t  Jo 


log(s) 


A:  = 0,  r < s 


^(r)l*l(l-s2W)  t^o,r< 


Substitute  this  result  into  equation  2.62  and  perform  algebraic  simplifications  using 
the  properties  of  Bessel  functions  expressed  in  equations  2. 7, 2. 8,  2.9  and  2.5  to  obtain 
the  result.  ■ 

Proposition  22  The  coefficients  Giik,m,n  are  given  by 


G Hk,m,n  — 


{^n,\k\)‘ 


I 1 

Jo  i-  — S 


(2.67) 
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Proof:  Recall  from  Proposition  14  in  the  discussion  of  the  Reisz  Decomposition  that 
can  be  expressed  as  the  integral 


2 

G'/ip,o(r,'0)  = - / , G{r  cos 'ip,r  sin  Ip;  p,t) 

7T  1 - _ ^2 


\/l  - 


dt. 


(2.68) 


Using  this  integral  representation,  we  can  invoke  Theorem  8 to  see  that  the  series 
coefficients  Gpk,m,n  are 


{•2-k  pi  pi  r\/l 

Jo  JoJ-V-s/l 


7T 


^ G{r  cosip,r  sinip;  p,t) 


1 — 


Tm{p)rJ\k\{rZn,ik\)e  dtdpdrdip. 

(2.69) 


We  can  express  the  iterated  integral  over  p and  t as  an  area  integral  over  the  disk  by 
introducing  polar  coordinates, 

p = scos{4>) 
t = ssin(0). 

Under  this  change  of  variable  the  coefficient  Gpk,m,n  can  be  written  as 

r>27r  pi  p2-K  /•! 


Y^Tm{scoscP)J\ki{rZn,\k\)e 
G{r  cos  ip,  r sin  ip;  s cos  (p,  s sin  cp)  ds  dcp  dr  dp;- 


—ik4> 


The  circular  symmetry  of  G can  be  used  to  rewrite  this  as 

ry  p2-K  pi  p2l\  pi 

ZCk,m,n 


III  / 1 2 

Jo  Jo  Jo  Jo  i-  — S 
G[r  cos{ip  — (p),r  sin(V’  - (p);s,  0)  ds  dcp  dr  dip. 


—ikip 


The  integral  can  be  grouped  as  follows, 
2ct. 


7T 


bV  rr,  , ..  f'P"  , , 

7 / y ^lm[s  cos  (p)  / / rJ\k\{rZn^lk\) 

Jo  Jo  ^ ^ Jo  Jo 


G{r  cos{ip  — (p),r  sin('0  — cp);s,  0)e  dip  dr  dcp  ds. 


A simple  change  of  variable  yields 

‘^C-k^rn, 


7T 


p\  p27T 

s cos  cp)e~^^^  dcP  j / rj\k\{rzn,\k\) 
Jo  Jo 


G{r  cos  Ip,  rsinip;s,Q)e  dip  dr  ds. 
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Combine  Proposition  21  with  the  definition  of  Tm^k{s)  in  equation  2.50  to  obtain  the 
proposition.  ■ 

As  the  following  corollary  shows,  Proposition  22  allows  the  structure  of  the 
matrix  Gt . . to  be  characterized. 

A,,  , 

Corollary  23  The  following  coefficients  are  zero, 

To  \k\  > m 

Gnk,m,n  = < 0 /c  even,m  odd  (2.70) 

I 0 /c  odd,  m even 


Proof:  Proposition  22  shows  that  Gfik,m,n  = 0 if  Tk^rn  = 0.  The  k,  m for  which  this 
occurs  are  identified  in  Proposition  18.  ■ 

We  will  now  compute  the  single  integral  in  equation  2.67.  Define  a function 
qk,m{s)  by 

qk,m{s)  = Y^Tk,m{s).  (2.71) 

It  follows  from  Proposition  19  that  for  m ^ k,  qk,m{s)  is  a polynomial.  Thus, 
equation  2.63  shows  that  the  coefficients  G^k,m,n,'m  k depend  on  the  Bessel  ex- 
pansion coefficients  of  polynomials.  The  following  proposition  states  a fact  from 
Bowman  [4]  that  leads  to  an  explicit  formula  for  such  expansions. 


Proposition  24  The  following  identity  holds  for  j,  k > 0, 


[ s^+'+"Vfe(zs)  ds  = V(-2) 

•^0  i=0 


— 1: 


U - i) 


•\1 


(2.72) 


This  proposition  is  immediately  applicable  to  the  computation  of  the  coeffi- 
cients G fXk,m,ri' 

Corollary  25  Let  m > k.  Let  Bij  be  the  lower  triangular  matrix  where 


= {-2y 


(j  - 0-' 


if  i < j, 


(2.73) 
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let  j{n,k)  be  the  vector  defined  by  ji{n,k)  = :!wi‘+i+djnM)  ^ q(m,k)  be  the 

vector  such  that  qfim^k)  is  the  coefficient  0/ in  the  polynomial  qk,m{s)-  The 
series  coefficients  of  the  Green’s  potential  are 


Gpk,m,n  = j(n,  kffiBq{m,k). 


(2.74) 


Proof:  Write  the  polynomial  qm,k{s)  as  + a2s'^^^  + + . . . + 0^-15'"“^ 

This  gives  the  elements  of  the  vector  qm,k  = [«i  «2  • • Apply  proposition  24  to 
obtain 


[ qm,k{s)Jk{zs)ds  = Y]aj  [ Jk{zs)  ds 

Jo  ■ Jo 


y-oVr-h+i) — I 


j i=0 


j i=0 

= j{n,k)'^Bq{m,k). 


(j  - 0! 


rj 


CHAPTER  3 

ALGORITHM  DEVELOPMENT 
3.1  Introduction 

In  this  chapter  we  develop  Proposition  9 into  an  algorithm  to  recover  the  den- 
sity function  from  scanner  data.  The  output  will  consist  of  (1)  the  series  coefficients 
of  the  recovered  density  function  and  (2)  the  recovered  density  function  sampled  on 
a radial  grid.  The  remaining  sections  of  this  chapter  detail  the  four  stages  of  the 
general  algorithm: 

• Computation  of  matrix  coefficients. 

• Data  Transformations. 

• Solution  of  the  system  of  equations. 

• Visual  reconstruction  on  a grid. 

Throughout  these  sections  it  is  assumed  that  routines  are  available  to  evaluate  Bessel 
functions  and  to  implement  the  East  Eourier  Transform  (EET),  such  as  those  de- 
scribed in  Press  et  al.  [31]. 

3.2  Matrix  Coefficients 

Proposition  9 shows  that  the  series  expansion  of  the  density  function  is  related 
to  the  series  expansion  of  the  data  by  a matrix  equation  where  the  elements  of  the 
matrix  are  J\^j^i{zn,\k\)P-k,m,n-  This  can  be  written  using  the  Riesz  Decomposition 

).  This  section  addresses  the  numerical  computation 
of  G-k,m,n  and  H_k,m,n,  truncated  at  m < M,\k\  < K and  n < N^. 
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The  series  coefficients  H^  rn^n  can  be  quickly  and  accurately  computed  during 
execution  using  the  formula  in  Corollary  17. 

Although  Proposition  24  presents  a fast  and  exact  formula  for  GiXk,m,ni  it  is  not 
well-suited  for  computational  purposes.  Term-by-term  integration  of  the  polynomial 
Qm,k{s)  relies  on  careful  cancelation  among  very  large  terms.  A slower  but  more  stable 
method  is  to  numerically  integrate  the  single  integral  in  Proposition  22.  Since  it  is 
expected  that  the  algorithm  will  be  used  to  recover  many  different  density  functions, 
it  is  efficient  to  precompute  and  store  these  coefficients. 

Equation  2.67  shows  that  involves  the  computation  of  {M  -|-  l)A^fc  inte- 

grals of  the  form 


Let  I be  even,  A = y and  sj  = jA,  for  j = 0, . . . , 1.  The  composite  Simpson’s  rule 
approximates  Im,n  by 


(3.1) 


where 


(3.2) 

(3.3) 


(3.4) 


(3.5) 


This  approximation  can  be  rewritten  in  vector  notation, 


= l/‘f'(0)  4/W(si)  2/W(sc)  ...  4/W(s,_.)  /W(l)l 
= fa'*>(0)  »<*'(*.)  . . . 9<‘>(si-i)  <ti‘'(l)l^ 

/(fc)  ~ fW  Jk) 

-‘m,n  J m ju  * 


(3.6) 

(3.7) 


(3.8) 
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Thus,  the  matrix  of  integrals  associated  with  the  trigonometric  coefficient  is 


Jo 

= : 


l(k)  ^ where 

r /■(*)  1 


(3.9) 


(3.10) 


L J M . 


(3.11) 


This  approach  exploits  the  form  of  the  integral  in  equation  3.1  by  minimizing  the 
number  of  function  evaluations.  A naive  application  of  Simpson’s  rule  to  each  integral 


real  impact  since  these  are  fairly  intensive  function  calls.  In  addition  to  the  function 
calls,  the  algorithm  involves  0((M  + l)NKl)  multiplications. 

Observe  that  requires  that  Tm,k{sj)  be  evaluated  for  k = 0, K,m  = 
0, . . . , M.  This  can  be  done  using  the  recursion  in  Proposition  20,  but  care  should 
be  taken  to  do  it  efficiently.  For  example,  Ti^i(sj)  should  be  stored  so  that  it  does 
not  have  to  be  recomputed  to  find  Ti^3(sj).  On  way  to  do  this  can  be  explained  by 
introducing  a triangular  matrix  T with  element  Tm^k  equal  to  T^,fc(s).  It  is  easily 
verified  that  the  polynomials  can  be  evaluated  without  redundancy  by  sequentially 
computing  the  diagonals  of  T.  Once  an  efficient  means  of  computing  Tm,k{s)  is 
identified,  we  can  find  the  matrix  elements  using  the  integration  method  described 
above. 


would  require  1{M  + 1)N  evaluations  of  the  form  fm\s)g^\s)  rather  than  only 
1{M  + 1)  evaluations  of  fm\s)  and  IN  evaluations  of  the  type  gn\s).  This  has  a 


3.3  Data  Transformations 


This  section  examines  the  transformations  necessary  to  compute  bk^rn  from 
scanner  data.  The  approach  described  below  can  be  illustrated  by  the  following 
path:  Tube  Data  — > Arc  Data  Trig.  Coefficients  — > Trig.-Tchebyshev  Coefficients. 
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3.3.1  Transforming  Tube  Data  into  Arc  Data 

Here  we  discuss  the  details  of  implementing  the  linear  transformation  M that 
was  introduced  in  section  1.3.  Recall  from  the  chapter  1 that  the  matrix  that  maps 
tubes  to  arcs  has  on  the  order  of  D'^  elements.  This  prohibitive  size  means  that  it  is 
preferable  not  to  explicitly  form  the  transformation  matrix  M defined  in  equation  1.9. 
Instead,  we  will  compute  the  arc  data  recursively.  Let  us  assume  that  the  tube  counts 
are  stored  in  a symmetric  D x D matrix  defined  by 

Uij  — mean  detections  in  tube  didj,  i,  j = 0, . . . , D — 1.  (3.12) 

The  arc  data  will  be  stored  in  the  D x D matrix  defined  by 

A,j  = mean  detections  in  arc  djdp+i)  mod  d • • • mod  /?,  i,  j = 0, . . . , D - 1. 

(3,13) 

The  following  proposition  provides  a means  of  transforming  tube  data  to  arc  data 
without  explicitly  multiplying  by  M. 

Proposition  26  The  elements  of  A satisfy  the  recursion, 

D 

^i,0  — Ud,i,  i = 0...D—l 

d=0 

i+j 

T -^i+j+1  mod  D,0  ^ ^ Ud  mod  D,i+j+l  mod  Di  J ~ 0,  . . . , D 2. 

d=i 

Proof:  The  formula  for  Ai^  follows  directly  from  the  definition  of  M.  j > 0 

refers  to  the  mean  number  of  emissions  such  that  at  least  one  photon  is  detected  by 
one  of  di,  c?(j+i)modo> ' ' • ’ ^(*+1+1).  mod  d-  This  can  be  partitioned  into  the  sum  of: 

• the  mean  number  of  emissions  such  that  at  least  one  photon  is  detected  by  one 

of  di,  mod  D 

• the  addition  detections  added  to  the  mean  by  including  detector  d(^i+j+i)^  mod  d- 
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The  first  summand  is  the  definition  of  Aij  while  the  second  summand  is  precisely 

■^i+j+l  mod  D,0  ES  u,  mod  D,t+j  + l mod  D-  ® 

There  are  two  points  that  can  significantly  affect  the  performance  of  an  im- 
plementation based  on  Proposition  26: 

1.  Explicit  use  of  the  modulo  operator  can  be  avoided. 

2.  E2  V,  mod  D can  be  computed  efficiently  from  previously  computed  par- 
tial sums. 


To  make  the  second  point  clear,  define 

i+j 

mod  Dji+j+l  mod  O' 

d=i 

Proposition  26  can  be  rewritten  in  terms  of  a recursion  for  both  A and  W as  follows, 


Corollary  27  The  first  column  of  A and  W are  given  by 

D 

Aifi  — Ud,ii  i = 0 . . . D — 1 

d=0 


IE, 


,,o 


i <D  ~l 
Uq,d-i  i = D — 1 


(3.15) 

(3.16) 


The  remaining  columns  (j  = 1, . . . , D — 1),  can  be  computed  recursively  by 


W = I ^*+15-1  + i < D - 1 

I Unv-i  + ^ = D-l 

= ^i,j  + Ai+j+\  mod  D,0  — j = 0,  . . . , D — 2. 


(3.17) 

(3.18) 


This  corollary  indicates  that  when  computing  row  j of  A,  the  column  of  W (rather 
than  the  entire  matrix)  should  be  stored  in  a vector  so  that  it  can  be  used  for  the 
[j  -I- 1)®^  column  of  A.  It  is  evident  that  such  an  implementation  would  involve  only 
0{D^)  additions. 
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Now  that  we  have  a means  of  computing  the  arc  data,  it  will  be  useful  to 
relate  the  elements  of  A to  the  arc  parameterization  described  in  section  1.3  (see 
Figure  1.2). 


Proposition  28  Aij  = b{pj,  9ij)  where  pj  = cos  ^{j  + 1)  and  6ij  = ^{j  + 1)  + 
l,j  = 0,...,D-1. 


Proof:  Since  each  detector  is  obtained  from  equi-spaced  partitions  of  the  unit  circle 
into  D subarcs,  detector  di  is  identified  with  the  arc. 


Aij  refers  to  the  mean  number  of  emissions  such  that  at  least  one  photon  is  detected 
by  one  of  di,  d(^i^i)  d,  ■■■,  mod  d-  These  detectors  form  an  arc  on  the  unit 
circle  consisting  of  the  points 


It  is  easily  verified  that  an  arc  with  endpoints  /?  > a is  parameterized  by 

p = cos  and  9 = Set  a = and  /3  = + j + 1)  to  obtain  the  theorem. 

■ 

3.3.2  Obtaining  Series  Coefficients 

Proposition  28  shows  that  column  j of  matrix  A contains  equally  spaced  an- 
gular samples  of  b{p,9),  shifted  by  ^j.  We  will  see  that  the  rows  of  A also  contain 
equally  spaced  samples  under  an  appropriate  change  of  variable.  This  sampling  pat- 
tern can  be  exploited  to  obtain  fast  approximations  to  the  series  coefficients  for  the 
data.  The  notation  in  the  following  proposition  will  be  simplified  with  the  introduc- 
tion of  the  Kronecker  delta. 


(3.19) 
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Proposition  29  A trapezoid  rule  approximation  to  the  trigonometric-Tchebyshev  co- 
efficients, bk^mt  Ihe  series  expansion  of  b{p,  6)  leads  to  the  double  summation 


bk,m  — 


2-6, 

D2 


mO 


D-1  D-\ 
j=0  1^0 


7T 


m-0  + 1 


D-e 


where 


\ j = D-l 

\ Aij,  otherwise 


(3.20) 


(3.21) 


Proof:  Theorem  6 tells  us  that  the  trigonometric-Tchebyshev  coefficients  of  b{p,  6) 
are  given  by 

bk,m  = Ck,m  f r b{p,  de  dp.  (3.22) 

J-l  Jo  y/l  — p2 

Combining  the  definition  of  Tm{p)  from  equation  2.1  with  the  change  of  variable 
(j)  = cos“^  p allows  us  to  express  equation  3.22  as 


^7T  ^27T 

/ / b{cos(j),9)  cos{m(f))( 

Jo  Jo 


b 


k,m  Ck^m 


,—ik9 


dO  d(p. 


(3.23) 


Apply  the  change  of  variable  9 = 9 — (j)  and  use  the  periodicity  of  the  integrand  to 
yield, 

pTT  /*2tT 

bk,m  = Ck,m  / b{cos(!),9 d9  d(j).  (3.24) 

Jo  Jo 

Approximate  the  inner  integral  by  means  of  the  composite  trapezoid  rule  to  obtain, 

D-l 


C/c,; 


27T 


/»7r  ^ ^ 

/ XI  Kcos  (j), 
1=0 


2n 


/ -I- (^)  cos(m0)e  d(f). 


(3.25) 


Note  that  equation  3.25  used  the  fact  that  the  integrand  is  periodic  with  period 
27t.  The  composite  trapezoid  rule  can  be  applied  again  to  the  remaining  integral  to 
obtain. 


1 1 


^/  + ^j)cos(m^i)e  (3.26) 
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The  terms  involving  ^ simply  give  the  first  and  last  terms  in  the  trapezoid  rule  the 
proper  factor  of  However,  equation  3.26  can  be  simplified  since  b{l,9)  = 0.  Thus, 
the  first  term  of  the  trapezoid  rule  may  be  omitted  since  it  is  0, 


The  proposition  follows  from  the  definition  of  Ck,m  in  Theorem  6.  ■ 

We  quote  a result  in  Stoer  [34]  that  will  allow  us  find  the  error  in  the  approx- 
imation used  in  Proposition  29, 

Theorem  30  Let  f € C^[a,  b],  h = and  Xj  = a + jh  for  j = 0,1, . . . ,D.  There 
exists  a £ {a,  b)  for  which  the  error  in  the  composite  trapezoid  rule  can  be  written, 


A standard  extension  of  Theorem  30  to  the  case  of  double  integrals  is  needed  so  that 
it  can  be  applied  to  trigonometric-Tchebyshev  expansions. 

Proposition  31  Let  / : 3?^  — >■  3?  6e  a function  with  continuous  second  partial  deriva- 
tives. Let  h = and  define  Xi  = a hi, i = 0,1, ...  ,n.  Similarly,  let  k = and 

define  = c + kj,  i = 0,1, . . . ,m.  There  exists  {r]i,p,i),  {r)2,  p.2)  £ [a,  b]  x [c,  d]  such 
that  the  composite  trapezoid  rule  applied  to  the  partition  {xi,yj}  has  error  term. 


Finally,  invoke  Proposition  28  to  arrive  at. 


(3.28) 


We  can  apply  Proposition  31  to  the  double  integral  approximation  in  Proposition  29 
to  obtain  the  following  result. 
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Corollary  32  The  error  in  the  trigonometric-Tchebyshev  coefficients  approximated 
in  Proposition  29  can  he  expressed  as 

tt'^ 

^k,m  ^fc,m  ~ ('k,m  qjj2 

where 

/ o2  L OL  \ 

a = - 2ik~{pi,0i)  ~ k‘^b{pi,ei)j 

pPh 

b = ( - — (p2,02)sin(cos“^(p2)) 

+^(p2,6'2)(2msin(cos~^(p2))sin(mcos“^(p2))  - pcos(p2)) 
-h{p2,92)m^Tm{p2)  ) 

for  some  (pi,  0i),  (p2,  ^2)  e [-1,1]  x [0,27rj. 

It  is  evident  that  the  method  of  Proposition  29  can  be  a poor  approximation 
for  large  m or  k.  We  outline  a method  described  in  Stoer  [34]  to  force  the  approxi- 
mated coefficients  to  exhibit  the  correct  asymptotic  behaviour.  The  approach  is  to 
approximate  the  true  coefficients  by  the  coefficients  of  an  interpolating  function.  The 
following  theorem  is  adapted  from  Gautschi  [12], 

Theorem  33  Let  T be  the  set  of  absolutely  continuous  21^ -periodic  functions  f : 
^ 5R.  Let  S : 3?^  be  a periodic  cubic  spline  interpolant.  Define  xi  = 

There  exist  attenuation  factors  Tk{N),k  G Z such  that  the  k^^  Fourier  coefficient  of 
the  approximating  function,  S {[f{xo)f{xi)  ■ ■ • f{xN-i)]') , is  f{xi)e-^'^^‘ . 

This  result  has  an  immediate  extension  to  two  dimensions. 

Corollary  34  Let  T2  be  the  set  of  2t: -periodic  and  absolutely  continuous  functions 
/ : 3?^  ^ R.  Define  Xi  = ^1,1  = 0, . . . , N - 1 and  yi  = ^1,1  = 0, . . . , M - 1.  Let 
S : T2  he  a periodic  cubic  tensor  spline  interpolant.  Denote  by  s{x,y)  the 

function, 


/ 

f{xo,yo)  • 

• f{xo,yM) 

\ 

V 

_ f{^N,yo)  • 

• f{^N,yM) 

/ 

S 


43 


that  approximates  / G based  on  the  sample  points  {xn,ym)i'>T'  = 0,  ■ ■ ■ , N,m  — 
0, . . . , M.  There  exist  attenuation  factors  Tk{-),  k E Z such  that 


s{x, 


dxdy  = Tm{M)Tk{N) 


1 


N-l  M-1 


j=0  1=0 


This  two  dimensional  corollary  is  relevant  to  the  computation  of  trigonometric- 
Tchebyshev  coefficients.  However,  its  application  is  not  as  direct  as  it  may  appear  at 
first  glance.  It  is  easy  to  show  that  the  trigonometric-Tchebyshev  coefficients  of  b{p,  9) 
are  equal  to  the  double  trigonometric  coefficients  of  b{cos(l),9),  which  the  corollary 
can  be  used  to  estimate.  However,  the  terms  of  the  2-dimensional  discrete  Fourier 
transform  of  6(cos  </),  9)  are  not  available.  The  geometry  of  the  scanner  restricts  us  to 
sample  points  that  were  calculated  in  Proposition  28.  These  samples  involve  shifts 
in  9 that  prevent  us  from  using  this  direct  approach.  Instead,  we  can  define  new 
functions  for  which 


1.  the  DFT  sample  points  are  available  and 

2.  the  double  trigonometric  coefficients  are  related  to  the  trigonometric-Tchebyshev 
coefficients  of  b{p,  9). 

The  following  theorem  takes  this  approach. 


Proposition  35  Let  the  functions  be  defined  by 


a\  — I b{cos  (j),9  + 4>)e  0 G [riTr,  nvr  -f  tt],  n even 

^ ' \ 6(cos((/)  — 7t),  0 -I- (^  — 0 G [nyr,  nrr -h  tt],  n odd 

(3.29) 

The  trigonometric-Tchebyshev  coefficients  ofb{p,9)  are  related  to  the  double  trigono- 
metric coefficients  of  9)  by 


U _ i 9 k)n  ^ 0 

1 145  ™ = 0 


(3.30) 
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Let  S : — )•  T2  denote  periodic  cubic  tensor  spline  interpolation.  Denote 

by  G the  matrix  with  elements  defined  by  = g^'^\^l,^j),l  = 

0, . . . , 2D  — 1,  j = 0, . . . , D — 1.  Let  bk^m  denote  the  coefficient 

1 p2TT  p2lT 

J j {SM^^\(f),e)e-^^^e-^'^’>>d9dp. 

There  exist  attenuation  factors  Tj{2D),Tj{D)  such  that  bk,m  = Tm{2D)Tk{D)bk,m- 

Proof:  Note  that  9)  = g{(f),  9).  Recall  that  the  formula  for  the  trigonometric- 

Tchebyshev  coefficients  is 


/I  p27T 

.1 


^k,m  — ; 


^-ike  Trnjp) 


Apply  the  change  of  variable  (j)  = cos  ^ p to  obtain 

C-K  /•27T 


^k,m  — ('k,r, 


b{cos(f),9)e  cos{m(j))  d9  d(f). 


0 JO 


(3.31) 


(3.32) 


A second  change  of  variable  9 = 9 — fi  yields 

"TT  f2-K 


bk  • 


^k,m 

('kfTn 

Ck 


b{cosfi,9  + (j))e  cos{m(f>)  d9  dfi 


'0  Jo 

r-n  p2-K 


g^^\(j),9)e  cos{m(j))  d9  d(j) 

p'Z77  P'ZTT 

r d9d(l). 

- Jo  Jo 


'0  Jo 

'>27t  />27r 


(3.33) 


The  stated  relationship  between  the  trigonometric-Tchebyshev  coefficients  of  b(p,  9) 
and  the  double  trigonometric  coefficients  of  g(fi,  9)  is  now  clear.  Apply  Corollary  34 


..  0-12D-1 

9l^l  = rk(n)U2D)^;^J^g('^)( 


to  see  that  there  exists  attenuation  factors  such  that 

a“’( Jpi,  -()e-~6>e-‘T'.  (3.34) 

1=0  j=0 

The  inner  sum  can  be  written  in  terms  of  Aij  to  obtain 

..  0-1  /D-l 

9f}n  = ^k{D)Tm{2D)^  ( Y cos(mTr) 


D-l 

+ 

1=1 


1=0  \j=l 


rkiD)Tm{2D)^Y^Aje  ''=S^'cos(m^j)e 


(=0  1=1 


(3.35) 
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It  follows  immediately  that  = Tk{D)Tm{‘2D)bk,rn-  ■ 

In  other  words,  we  have  bk,m  = gi%  « Tk{D)Tm{2D)bk,m  = h,m-  This  is  the 
method  that  will  be  used  to  approximate  the  series  coefficients  of  the  data  in  the 
remainder  of  this  work.  The  examples  will  use  periodic  cubic  spline  approximations 
for  which  the  attenuation  factors  are  calculated  in  Stoer  [34], 


rj{N)  = 


sm 


ni 

N 


^ / l + 2cos2^ 


(3.36) 


We  proceed  to  discuss  an  efficient  implementation  of  this  approach.  bk,m  can  be 
computed  quickly  since  the  summation  over  I in  equation  31  is  precisely  the  discrete 
Fourier  transform  and  the  summation  over  j is  a shifted  version  of  the  discrete  cosine 
transform.  It  is  thus  immediately  evident  that  the  inner  summation  in  the  formula 
for  bk^m  can  be  implemented  as. 


bk,m  = 5]  COS  (m^{j  + 1))  (3.37) 

j=0 

where  FFT(74)fcj  denotes  the  matrix  whose  column  is  the  fast  Fourier  transform  of 
the  column  of  A.  The  next  proposition  will  aid  us  in  devising  a fast  implementation 
for  the  remaining  summation. 

Proposition  36  Let  Vj,j  = be  a vector  of  length  D.  Let  v be  an  even 

extension  of  v defined  by 


D-l 


^3  = 


0 

i = 0 

^3 

j = 1 

2vd 

j=.D 

. ^3-D 

j = D + 1,...,2D-1 

'm  of  V 

is  equal  to 

^^cos  (m^(j  + l))u,- 

3=0 

Proof:  The  element  of  the  discrete  Fourier  transform  of  v is 

2D-1 

2D 


^ 2D-1 

FFT(5)„  = 


3=0 


46 


-D-l  2D-1 

y InJ 


Substituting  the  definition  of  v into  this  equation  gives 

= ^(E 

\j=o 

This  simplifies  to 

1 / . . . . 
FFT(fi)^  = — 2vd  cos(mTr)  + ^ 


Wj-oe  -|_  2uj5  cos(mTr) 

j=D+\ 


j=o 


This  proposition  can  be  used  to  compute  equation  3.37  if  we  assign  vj  = 
j = 1, . . . , J9  — 1.  A MATLAB  algorithm  to  compute  the  co- 
efficient approximations,  Tm{2D)Tk{D)hk^rn  using  fast  Fourier  transforms  is  in  the 
apendix  (function  coeff).  Its  computational  complexity  is  0{D'^  log  D). 


3.4  Solution  of  the  System  of  Equations 


The  density  coefficients  satisfy  the  infinite  dimensional  matrix  equations  given 
in  Proposition  9.  It  is  assumed  that  we  can  identify  K and  Nk,  M^,  k = —K, . . . ,K 
such  that  we  can  approximate  the  trigonometric-Bessel  coefficients  of  A(p,  6)  by  solv- 
ing the  truncated  systems 

^ l(k)^  ^3_3g^ 


where 

lik) 

and 


bk,o 

, = 

1 

1 

^k,Nk 

(3.39) 


p(fc) 


Jfk\  + li^l,\k\)P-kfi,l 


'^\k\  + li^Nk,\k\)P-k,0,Nk 


IT 


(3.40) 


_ J\k\+li^l,\k\)P-k,MkA  •••  J^k\  + li^Nk,\k\)P-k,Mk,Nk  _ 

It  will  be  convenient  to  write  these  2K  + 1 matrix  equations  as  a single  block  diagonal 
system  of  equations, 


PX  = b 


(3.41) 
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where 


- li-K)  - 

■ A(-^)  ■ 

- p{-K) 

0 

b = 

,A  = 

, and  P = 

l(K)  _ 

A(^) 

0 

p{K) 

It  can  be  verified  numerically  that  the  singular  values  of  P decay  very  rapidly  to 

zero.  Consequently,  we  need  to  provide  additional  information  regarding  A in  order 

to  find  a stable,  unique  solution.  Specifically,  we  will  determine  decay  rates  for  the 

coefficients  of  A under  reasonable  smoothness  conditions  and  then  use  a regularization 

method  to  encourage  the  solution  to  decay  similarly. 

3.4.1  Decay  Rates 

We  invoke  the  following  result  from  Watson  [40]  to  determine  the  rate  of  decay 
of  the  Bessel  coefficients. 


Theorem  37  Suppose  that  the  density  X{r,6)  is  such  that  the  functions  A*,(r)  = 
\/r  X{r,9)e~'^^^  dd  are  of  bounded  total  variation.  Then  the  coefficients  in  the 

order-k  Bessel  expansion  of  the  function  Xk{r)  satisfy  Xk,n  < 

^n,k 

The  requirements  of  Theorem  37  seem  to  imply  that  the  density  should  approach  0 
at  the  origin.  The  next  proposition  shows  that  we  can  weaken  this  assumption  to  the 
more  appropriate  restriction  that  the  density  is  constant  on  a tiny  disk  centered  at 
the  origin. 


Proposition  38  If  the  density  is  a disk  about  the  origin  defined  by 

Hp  { 0,  otherwise  ’ 
then  its  trigonometric-Bessel  coefficients  are 

X{r,^)e-^^‘^Jk{rZn,\k\)rdrd(f)=  I | ^ ^ ^ n ' 

Proof:  For  A:  0 the  result  is  clear.  For  A;  = 0 use  equation  2.7  to  perform  the 

integration  to  obtain  The  result  is  a consequence  of  the  approximation 

\Jp{zn,pr)\  < ^^0=  from  Bowman  [4].  ■ 
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Theorem  37  and  Proposition  38  motivate  us  to  regularize  the  solution  of  3.41 
to  encourage  the  solution  to  exhibit  these  decay  rates.  An  appropriate  regularizing 
operator  is  therefore  defined  by 


r 


L = 


0 


0 


(3.43) 


where 


0 


. 0 izN,,\k\)^ 

We  can  transform  3.41  into  standard  form  using  P = 
problem  becomes  one  of  regularizing 


(3,44) 

PL~^,\  = LX  so  that  the 


PA  = b.  (3.45) 

3.4.2  Regularization 

This  section  addresses  the  regularized  solution  phase  of  the  reconstruction  al- 
gorithm. There  are  many  reasonable  approaches  to  computing  a regularized  solution 
to  equation  3.45  - Tikhonov’s  method,  truncated  SVD,  etc..  We  choose  the  conjugate 
gradient  method  because  of  its  speed,  ease  of  parameter  selection,  and  most  impor- 
tantly its  good  performance  on  test  densities.  These  issues  will  be  addressed  later  in 
this  section  and  in  chapter  4.  Unfortunately,  the  regularizing  properties  of  the  CGM 
are  not  as  well  understood  as  the  classical  approaches. 


Conjugate  Gradient  Method  The  conjugate  gradient  method  (CGM)  is  an 
iterative  optimization  algorithm  for  solving  the  quadratic  problem, 

minimize  f{x)  = ^x^Ax  - y^x,  (3.46) 

where  A is  positive  definite.  A derivation  of  the  CGM  and  a discussion  of  its  conver- 
gence properties  can  be  found  in  Bertsekas  [2]  and  Luenberger  [25].  The  CG  iteration 
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is 


Pk 

<r- 

Adk, 

^fc+1 

<r- 

Xk  + Cikdk,  where  > 

d>kPk 

Tk+l 

Tk  - OLkPk 

dk+i 

t- 

Tk+i  + /dfcdfc,  where  (5k  = 

and  the  initial  values  are  xq  = 0,  do  = ’’o  = b. 

It  has  been  observed  by  Plato  [30],  Van  Der  Sluis  [37],  and  Vogel  [39]  that  the 
low  frequency  components  of  the  solution  converge  faster  than  the  high  frequency 
components.  Therefore,  the  CGM  can  be  used  to  obtain  regularized  solutions  of  3.45 
by  applying  it  to  the  associated  normal  equations.  Specifically,  replace  A and  y 
in  equation  3.46  by  P^P  and  P^b,  respectively.  The  degree  of  regularization  is 
controlled  by  the  number  of  iterations. 

Hansen  [17]  shows  how  to  use  the  results  in  Golub  [14]  to  put  the  T'*  iterate 
of  the  CGM  into  an  optimization  framework  analogous  to  the  familiar  framework  for 
Tikhonov  regularization.  Specifically,  the  iterate  solves, 

min  \\PX  - 5||2  subject  to  A e )Ci{P^P,  P^b),  (3.47) 

where  K,i{P^P,  P'^b)  is  the  Krylov  subspace  span{P^6,  P^PP^b, . . . , 

Since  P is  block  diagonal,  equation  3.47  can  be  written 

K 

min  ^ _ 5^112  subject  to  A^^^  e K.i{P^'''>'^ P^^\  P^’^^'^b).  (3.48) 

k=-K 

It  is  thus  evident  that  the  iterate  of  the  CGM  applied  to  for  each 

k yields  the  same  solution  as  the  iterate  of  the  CGM  applied  to  PA  = b. 

It  is  an  interesting  observation  that  the  CGM  applied  to  P^PA  = P'^b  is  equiv- 
alent to  applying  the  preconditioned  conjugate  gradient  method  (see  Golub  [14])  to 
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P^PX  = P^b  with  preconditioner  L~^.  However,  the  normal  purpose  of  precondi- 
tioning is  to  accelerate  convergence,  whereas  our  motivation  was  to  regularize  LX. 

Bjorck  [3]  describes  how  to  apply  the  CGM  to  the  normal  equations  without 
explicitly  forming  the  matrix  product.  We  wish  to  apply  Bjorck’s  algorithm  to  prob- 
lem 3.45.  The  comments  above  indicate  that  this  can  be  achieved  by  applying  the 
algorithm  separately  to  each  of  the  2K  + 1 blocks.  An  efficient  implementation  must 
take  advantage  of  the  special  structure  of  these  blocks.  Recall  that  P^'^'>  = G + uv^ 
where  the  nonzero  rows  of  G are  identified  in  Corollary  23.  The  following  propo- 
sition shows  how  to  exploit  this  structure.  (The  order  of  operations  indicated  by 
the  brackets  is  crucial  for  an  efficient  implementation.)  Bjorck’s  algorithm  should  be 
modified  to  use  the  multiplication  outlined  in  Proposition  39. 

Proposition  39  Let  X € 3?^*  and  y G 3?^*.  Define  Mk  = ^ ^ 

sjfjMfcxAffc  satisfy  P = G + uv^ , where  Gij  = 0 for  i ^ k + 1,  k + 3, . . . , k + 2Mk  ~ 1- 
Define  a matrix  G that  contains  the  nonzero  rows  of  G by  Gtj  = Gk-i+2i,j,i  = 
1, . . . , Mk,j  = 1,  • • • 5 Nk-  Define  vectors  x and  y by 

yi  = yk-i+2u  fori  = l,...,Mk 

_ ^ i {Gx)i,  i = k-l  + 2j  for  j e {l,...,Mk} 

* [ 0,  otherwise 

The  multiplication  P^y  = G^y  -I-  v{u^y)  has  computational  complexity  0{MkNk  + 
2Nk)  and  the  multiplication  Px  = x + u{v^x)  has  complexity  0{MkNk  + 2A^*,). 
Parameter  Selection 

As  with  all  regularization  methods,  the  CGM  involves  a parameter  that  con- 
trols the  degree  of  smoothing.  In  this  case  the  parameter  is  i,  the  number  of  iterations. 
If  too  few  iterations  are  performed,  then  the  solution  will  be  oversmoothed.  On  the 
other  hand,  too  many  iterations  will  result  in  undersmoothing.  General  parameter 
selection  methods  are  described  by  Hansen  [18,  17].  Implementation  of  a black-box 
parameter  selection  method  is  problematic  (see  Hansen  [17]);  parameter  selection  is 


probably  best  done  with  some  human  guidance.  The  examples  in  chapter  4 will  use 
the  CGM  with  the  iteration  i chosen  by  visual  inspection  with  guidance  from  the 
L-curve  method.  Even  exculsive  use  of  visual  inspection  is  a feasible  approach  be- 
cause, in  theory,  at  most  maxiV^  iterates  need  to  be  examined.  More  importantly, 
the  examples  indicate  that  the  conjugate  gradient  method  performs  well. 


The  recovered  density  function  can  be  sampled  at  any  point  on  the  unit  disk 
by  means  of  formula  2.14.  However,  equation  2.14  is  computationally  expensive  since, 
for  every  sample  point,  we  must  accumulate  a double  summation  as  well  as  evaluate 
Bessel  functions.  This  involves  ~ 0{{2K  + l)iVo)  Bessel  evaluations 

and  multiplications.  Thus,  a naive  algorithm  to  sample  the  density  function  on  an 
L X L grid  would  have  {0){L'^{2K  -h  l)A^o)  multiplications  and  Bessel  evaluations. 
For  L large  enough  to  produce  a visually  pleasing  plot  of  the  density,  this  would  be 

the  most  time  consuming  phase  of  the  overall  algorithm. 

3.5.1  A Fast  Visual  Representation  Algorithm 

We  will  see  that  a significant  reduction  in  computational  complexity  can  be 
achieved  by  reconstructing  A(r,  6)  on  a radial  grid  with  equally  spaced  angular  sam- 
ples. Specifically,  a fast  algorithm  for  sampling  A(r,  ^j),j  = 0, 1, . . . , C = 2",  C > 
2K  will  be  derived. 

Let  / (r)  be  the  vector  of  Fourier  coefficients  determined  by 


3.5  Visual  Representation  of  an  Image 


fk{r) 


^n,k-NJ\k-N\{Zn,\k-N\'>')  N ~ K < k < N , (3.49) 


otherwise 


0<k<K 


where  k=0,. . . ,N-1.  Let  A(r)  be  a vector  of  angular  samples  at  radius  r defined  by. 
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This  leads  to  the  identity 

K Nk  N-l 

^j(0=  X]  V„J|fc|(z„,|fc|r)e*'=^”  = (3.51) 

k=—K  n=0  fc=0 

This  shows  that  the  vector  A(r)  = [Ao(r), . . . , Ac-i(r)]  is  precisely  the  discrete  inverse 
Fourier  transform  of  f{r).  Thus,  A(r)  can  be  efficiently  computed  using  the  inverse 
FFT, 

A(r)  = C IFFTiMr)).  (3.52) 

This  method  involves  Nk)  Bessel  function  evaluations  and  and  0(ClogC  + 

z2k^-K^k)  multiplications.  A reconstruction  on  the  unit  disk  can  be  obtained  by 
applying  the  algorithm  to  suitably  spaced  radii.  Furthermore,  if  many  density  func- 
tions are  to  be  reconstructed  on  the  same  grid,  the  Bessel  function  evaluations  need 
to  be  computed  only  once. 

3.6  A Reconstruction  and  Visualization  Algorithm 

The  results  of  the  preceding  sections  can  be  combined  to  arrive  at  the  following 
algorithm  for  numerical  reconstruction  and  visualization  of  the  reconstructed  density., 

• (Precompute  and  store  Gfik,m,n-) 

• Choose  K,Nk,Mk. 

• Compute  Hk,m,n- 

• Convert  the  tube  counts  to  arc  counts  using  function  xform. 

• Compute  the  trigonometric-Tchebyshev  coefficients  using  function  coeff. 

• Apply  the  transformation  P = PL~^,  where  P is  the  block  diagonal  matrix 
defined  by  equation  3.42. 

• Compute  a regularized  solution  to  the  transformed  block  diagonal  system  PA  = 
b using  function  cgm. 
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• Reconstruct  images  from  the  untransformed  iterates  using  function  recon. 

• Choose  the  best  image. 

We  can  combine  the  flop  counts  for  the  various  sections  of  the  algorithm  to  arrive  at 
the  total  computational  complexity.  If  we  denote  the  number  of  rings  composing  the 
reconstruction  grid  by  R,  then  the  total  complexity  of  computing  and  reconstructing 
all  I iterations  is 

K K 

0{DHogD  + IRC\ogC + IR  Nk  + I Y ^^Nk).  (3.53) 

k=-K  k^-K 

We  can  simplify  this  expression  a bit  if  we  assume  that  Nk  and  are  chosen  to  be 
constants, 

0{D^  log  D + IRC\ogC  + IRK  No  + IKMoNq).  (3.54) 

Furthermore,  I,  R,  K,  Nk,  Mk  are  all  typically  very  small  (e.g.  < 40). 

For  comparison,  the  comments  in  chapter  1 indicate  that  I iterations  of  the 
EM-ML  reconstruction  on  a grid  of  the  same  size  would  have  complexity 

0{IRCD^).  (3.55) 

It  should  be  noted  that  this  comparison  is  not  entirely  fair  since  a good  implementa- 
tion of  the  EM-ML  algorithm  can  take  advantage  of  the  sparsity  of  Q. 


CHAPTER  4 

NUMERICAL  EXPERIMENTS 
4.1  Introduction 

The  stated  goal  of  this  work  is  the  development  of  a fast  reconstruction  method 
that  incorporates  finite  detector  width.  The  only  other  approach  that  can  accommo- 
date finite  detector  width  is  the  discretization  approach  described  in  section  1.2.2; 
the  standard  solution  method  associated  with  this  approach  is  the  ML-EM  iteration. 
Eor  this  reason,  the  success  of  this  work  in  achieving  its  stated  goal  is  assessed  by 
comparing  reconstructions  using  the  orthogonal  function  approach  with  those  gener- 
ated by  the  ML-EM  algorithm.  This  evaluation  addresses  speed  and  image  quality 
and  is  based  on  a variety  of  test  density  functions,  called  phantoms.  Eor  ease  of 
exposition,  the  Truncated  Orthogonal  Series  reconstruction  method  described  herein 
will  be  referred  to  as  the  “TOS”  method. 

Eor  all  of  the  examples  that  follow,  it  is  assumed  that  the  scanner  ring  is 
comprised  of  128  detector  elements.  The  reconstructions  using  the  TOS  method 
set  K = 40,  N).  = 40  — \k\,k  = The  ML-EM  reconstructions  are 

performed  on  a standard  128  x 128  grid.  The  TOS  reconstructions  were  obtained 
using  the  MATLAB  algorithms  described  in  chapter  3.  The  ML-EM  reconstructions 
were  obtained  from  compiled  C code.  60  iterations  of  both  methods  were  computed, 
and  the  best  reconstruction  was  chosen  visually  for  inclusion  in  this  chapter. 

The  remaining  sections  show  that  the  image  quality  varied  among  the  different 
phantoms.  The  reconstruction  speed,  however,  is  not  a function  of  the  phantom  and 
so  we  discuss  it  now.  In  every  case,  the  time  it  took  to  perform  25  ML-EM  iterations 
was  within  a ±1  sec  of  528  sec  on  a SUN  Ultrasparc  workstation.  This  contrasts 
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Figure  4.1:  A smooth  radial  phantom  (cross-section  from  the  center  to  the  perimeter). 

dramatically  with  TOS,  which  always  took  under  4.8  sec  for  all  25  iteration.  Fur- 
thermore, these  results  are  significantly  biased  against  TOS  since  it  was  implemented 
using  a MATLAB  interpreter  rather  than  compiled  code. 


Phantom  Definitions  The  first  test  phantom  is  a radial  phantom  given  by 
the  equation 


A cross  section  of  this  test  phantom  that  begins  at  the  center  is  shown  in  the  first 
panel  of  figure  4.1.  The  constant  c was  chosen  so  that  a mean  of  two  million  emissions 
would  be  generated. 

Data  Generation  Following  the  Shepp- Vardi  model  expressed  by  equation  1.1, 
the  mean  number  of  emissions  in  tube  t can  be  calculated  by 
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4.2  A Smooth  Radial  Phantom 


(4.1) 
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(4.2) 
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The  mean  tube  counts  were  obtained  by  numerically  approximating  equation  4.2 
using  a standard  quadrature  rule.  The  advantage  of  experimenting  with  a radial 
phantom  lies  in  exploiting  the  circular  symmetries  shared  by  the  kernel  and  the  den- 
sity function.  For  a scanner  ring  with  128  detectors,  only  64  such  integrals  need  to 
be  computed.  The  remaining  mean  tube  counts  can  be  obtained  by  simple  rotational 
symmetries.  This  is  a feasible  way  to  generate  data  that  is  independent  of  the  recon- 
struction algorithm.  Once  the  mean  tube  counts  were  generated,  Poisson  noise  was 
added  to  obtain  a realistic  set  of  tube  counts. 

Results  Cross  sections  of  the  ML-EM  and  TOS  reconstructions  are  shown 
in  Figure  4.2;  the  first  column  contains  the  TOS  reconstructions  while  the  second  col- 
umn contains  the  ML-EM  reconstructions.  The  images  on  the  top  row  were  produced 
using  mean  tube  data.  Data  with  Poisson  noise  was  used  for  the  reconstructions  in 
the  second  row.  Both  methods  perform  quite  well  for  the  noise- free  data,  but  the 
ML-EM  method  exhibits  more  over-smoothing  for  radii  greater  then  0.6.  The  TOS 
reconstruction  degrades  a bit  near  the  perimeter  when  noise  is  added.  The  ML-EM 
reconstruction  degrades  significantly  when  noise  is  added.  It  can  be  stated  unequiv- 
ocally that  for  this  smooth  phantom  the  TOS  method  is  preferable,  particularly  in 
the  presence  of  noise.  Note  that  TOS  could  have  attained  better  results  if  we  had 
included  the  a priori  knowledge  that  the  phantom  was  radial.  This  could  have  been 
achieved  by  setting  /L  = 0.  Since  such  knowledge  is  not  normally  available,  a fair 
comparison  with  ML-EM  should  probably  not  make  use  of  it. 

4.3  Discontinuous  Phantoms 

Phantom  Definitions  In  this  section,  the  TOS  method  is  tested  on  the  3 
standard  phantoms  shown  in  figures  4.3  and  4.4.  The  simple  phantom  and  the 
Hoffman  phantom  have  been  used  extensively  in  the  literature  (e.g.  see  Mair  et 
al.  [26]).  All  of  these  phantoms  are  scaled  so  that  the  mean  number  of  emissions 
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x10®  (c)  X10^  (d) 


Figure  4.2:  Reconstruction  of  a smooth  radial  phantom,  (a)  Cross  section  of  a TOS 
reconstruction  using  noise-free  data,  (b)  Cross-section  of  an  ML-EM  reconstruction 
using  noise-free  data,  (c)  Cross  section  of  a TOS  reconstruction  using  noisy  data, 
(d)  Cross-section  of  an  ML-EM  reconstruction  using  noisy  data. 
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Figure  4.3:  Three  discontinuous  phantoms,  (a)  Simple  phantom,  (b)  Cardiac  phan- 
tom. (c)  Hoffman  phantom. 

is  two  million.  Since  these  phantoms  are  step  functions  on  a 128  x 128  grid,  they 
represent  a difficult  challenge  for  the  TOS  method.  On  the  other  hand,  this  fact 
makes  the  phantoms  particularly  well  designed  to  work  with  ML-EM.  Recall  from 
section  1.2.2  that  this  is  an  assumption  upon  which  the  ML-EM  method  is  based. 

Data  Generation  The  method  of  numerical  integration  is  no  longer  feasible 
for  general  nonradial  phantoms.  The  lack  of  symmetry  means  that  a double  integral 
has  to  be  approximated  for  all  128^  tubes.  Instead,  we  will  use  equation  1.4  to 
approximate  the  mean  tube  counts  since  Q and  A are  known.  The  data  that  we 
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Figure  4.4:  Cross-sections  through  the  center  of  the  phantoms,  (a)  Simple  phantom, 
(b)  Cardiac  phantom,  (c)  Hoffman  phantom. 
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will  take  to  be  noise  free  will  actually  have  errors  due  to  the  discretization  of  the 
kernel.  This  has  a two-fold  effect:  (1)  discretization  errors  associated  with  the  ML- 
EM  method  are  not  accounted  for  and  (2)  the  noise-free  TOS  reconstructions  are  in 
reality  not  quite  noise-free.  To  minimize  these  effects,  the  data  was  generated  on  a 
finer  grid  (256  x 256)  than  the  grid  that  was  used  for  reconstructions  (128  x 128).  This 
technique  is  convenient  and  consistent  with  our  philosophy  of  providing  a conservative 
appraisal  of  the  TOS  method. 

Results  The  reconstructions  for  the  Simple,  Cardiac  and  Hoffman  phantoms 
are  contained  in  figures  4. 5,4. 7,  4.9,  respectively.  Figures  4.6,  4.8  and  4.10  contain 
plots  of  horizontal  cross-sections  through  the  origin  for  each  reconstruction.  In  gen- 
eral, the  TOS  reconstructions  are  of  similar  but  distinctly  inferior  quality  than  the 
ML-EM  reconstructions.  That  is,  both  methods  were  able  to  identify  the  same  struc- 
tures, but  the  ML-EM  reconstructions  are  crisper  with  less  artifact.  Both  methods 
to  degrade  when  noise  is  added.  However,  the  nature  of  the  degradation  is  different; 
the  ML-EM  reconstructions  became  speckled  whereas  wavy  oscillations  showed  up  in 
the  TOS  images.  These  oscillations  were  particularly  evident  in  the  region  between 
the  phantom  and  the  scanner  ring.  The  line  plots  reveal  that  this  noise  is  substan- 
tially smaller  in  magnitude  than  the  phantom  itself,  and  so  it  could  be  successfully 
removed  with  some  simple  postprocessing.  Nevertheless,  it  should  be  noted  that  only 
for  the  cardiac  phantom  did  the  ML-EM  method  show  noise  in  the  region  between 
the  body  and  the  scanner  ring.  These  observation  should  not  come  as  a surprise  since 
the  ML-EM  algorithm  always  approximates  phantoms  by  stepwise  constant  functions 
whereas  the  TOS  method  is  not  as  well  suited  to  approximate  such  regions. 
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(a) 


(b) 


Figure  4.5:  Reconstruction  of  the  Simple  Phantom,  (a)  A TOS  reconstruction  us- 
ing noise  free,  (b)  An  ML-EM  reconstruction  using  noise-free  data,  (c)  A TOS 
reconstruction  using  noisy  data,  (d)  An  ML-EM  reconstruction  using  noisy  data. 
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Figure  4.6:  Cross-section  of  the  Simple  Phantom  reconstructions,  (a)  A TOS  recon- 
struction using  noise  free,  (b)  An  ML-EM  reconstruction  using  noise-free  data,  (c) 
A TOS  reconstruction  using  noisy  data,  (d)  An  ML-EM  reconstruction  using  noisy 
data. 
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Figure  4.7:  Reconstruction  of  the  Cardiac  Phantom,  (a)  A TOS  reconstruction 
using  noise  free,  (b)  An  ML-EM  reconstruction  using  noise-free  data,  (c)  A TOS 
reconstruction  using  noisy  data,  (d)  An  ML-EM  reconstruction  using  noisy  data. 
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Figure  4.8:  Cross-Section  of  the  Cardiac  Phantom  reconstructions,  (a)  A TOS  recon- 
struction using  noise  free,  (b)  An  ML-EM  reconstruction  using  noise-free  data,  (c) 
A TOS  reconstruction  using  noisy  data,  (d)  An  ML-EM  reconstruction  using  noisy 
data. 
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(a) 


(b) 


Figure  4.9:  Reconstruction  of  the  Hoffman  Phantom,  (a)  A TOS  reconstruction 
using  noise  free,  (b)  An  ML-EM  reconstruction  using  noise-free  data,  (c)  A TOS 
reconstruction  using  noisy  data,  (d)  An  ML-EM  reconstruction  using  noisy  data. 
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Figure  4.10:  Cross-section  of  the  Hoffman  Phantom  reconstructions,  (a)  A TOS 
reconstruction  using  noise  free,  (b)  An  ML-EM  reconstruction  using  noise-free  data, 
(c)  A TOS  reconstruction  using  noisy  data,  (d)  An  ML-EM  reconstruction  using 
noisy  data. 
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4.4  Conclusions 

The  examples  in  this  chapter  provide  a means  to  asses  the  potential  of  the 
orthogonal  function  method  of  reconstructing  PET  images  measured  against  the  cur- 
rent standard  ML-EM  algorithm.  As  discussed  above,  this  assessment  is  a conserva- 
tive one  since  the  non-radial  comparisons  are  biased  in  favor  of  the  ML-EM  method 
with  regard  to  implementation  platform,  data  generation  and  the  construction  of 
the  phantoms  themselves.  Nevertheless,  the  orthogonal  function  algorithm  produced 
promising  images  that  compared  favorably  with  the  ML-EM  images.  More  precisely, 
the  smooth  example  suggested  that  TOS  will  outperform  EM  on  noisy  phantoms  that 
are  well  represented  by  a truncated  trigonometric-Bessel  series  expansion  whereas  the 
discontinuous  phantoms  suggested  that  EM  will  outperform  TOS  on  phantoms  that 
are  difficult  to  approximate  with  a truncated  trigonometric-Bessel  series. 

ft  is  clear  that  the  TOS  algorithm  dramatically  outperforms  the  EM  algorithm 
in  terms  of  speed.  Although  acceleration  techniques  can  be  applied  to  the  ML-EM 
iteration  (see  Kaufman  [20]),  no  modification  of  the  ML-EM  method  has  approached 
the  speed  of  the  TOS  method.  In  fact,  it  is  unlikely  that  any  acceleration  method  can 
be  as  fast  as  TOS  since  the  above  examples  indicate  that  TOS  converges  in  less  time 
than  a single  EM-style  iteration.  A good  compiled  implementation  of  the  algorithm 
described  in  chapter  3 would  appear  nearly  instantaneous  to  the  end  user.  This  is  the 
only  reconstruction  method  that  accounts  for  finite-width  detectors  that  can  make 
such  a claim. 

In  short,  the  preferred  method  should  be  determined  by  the  nature  of  the 
density  function  and  the  importance  of  computation  time. 
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